Changeset 34404

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/psphot.h (modified) (3 diffs)
• src/psphotAddNoise.c (modified) (2 diffs)
• src/psphotApResid.c (modified) (1 diff)
• src/psphotBlendFit.c (modified) (7 diffs)
• src/psphotDeblendSatstars.c (modified) (3 diffs)
• src/psphotExtendedSourceAnalysis.c (modified) (1 diff)
• src/psphotExtendedSourceFits.c (modified) (1 diff)
• src/psphotFindDetections.c (modified) (2 diffs)
• src/psphotFitSourcesLinear.c (modified) (1 diff)
• src/psphotGuessModels.c (modified) (1 diff)
• src/psphotKronIterate.c (modified) (1 diff)
• src/psphotMagnitudes.c (modified) (1 diff)
• src/psphotRadialApertures.c (modified) (1 diff)
• src/psphotRadialProfileWings.c (modified) (2 diffs)
• src/psphotReadout.c (modified) (5 diffs)
• src/psphotSourceSize.c (modified) (1 diff)
• src/psphotSourceStats.c (modified) (3 diffs)
• src/psphotStackImageLoop.c (modified) (1 prop)
• src/psphotStackMatchPSFs.c (modified) (1 diff)
• src/psphotVisual.c (modified) (14 diffs)

trunk/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            psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources);
+bool            psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal, pmSourceMode showmode);
+bool            psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources, pmSourceMode showmode);
 bool            psphotVisualShowFlags (psArray *sources);
 bool            psphotVisualShowSourceSize (pmReadout *readout, 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);
@@ -511 +515 @@
     );
 
+bool psphotSatstarProfileModel (pmSource *source, psImageMaskType maskVal);
+bool psphotSatstarProfileCreate (pmSource *source, psVector **logRmodelOut, psVector **logFmodelOut, psVector *logR, psVector *flux, float Rmax);
+bool psphotSatstarProfileOp (pmSource *source, psImageMaskType maskVal, float FACTOR, pmModelOpMode mode, bool add);
+bool psphotVisualRadialProfileSatstar (pmSource *source, psImageMaskType maskVal);
+bool psphotAddOrSubSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex, psMetadata *recipe, bool add);
+bool psphotSatstarPhotometry (pmSource *source);
+
 bool psphotSourceMemory(pmConfig *config, const pmFPAview *view, const char *filerule);
 bool psphotSourceMemoryReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index);

trunk/psphot/src/psphotAddNoise.c

@@ -76 +76 @@
     }
 
+    psphotVisualShowImage (readout);
+
     // loop over all source
     for (int i = 0; i < sources->n; i++) {
         pmSource *source = sources->data[i];
+
+        // add or subtract noise for a saturated star.  satstars modeled as a radial profile
+        // need special handling
+        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) {
+            psphotSatstarProfileOp (source, maskVal, FACTOR, PM_MODEL_OP_NOISE, add);
+            continue;
+        }
 
         // skip sources which were not subtracted
@@ -92 +101 @@
     }
 
+    psphotVisualShowImage (readout);
+
     return true;
 }

trunk/psphot/src/psphotApResid.c

@@ -200 +200 @@
         if (source->mode &  PM_SOURCE_MODE_FAIL) SKIPSTAR ("FAIL STAR");
         if (source->mode &  PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR");
+
+        // skip saturated stars modeled with a radial profile 
+        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) SKIPSTAR ("SATSTAR");
 
         if (source->mode &  PM_SOURCE_MODE_EXT_LIMIT) SKIPSTAR ("EXTENDED");

trunk/psphot/src/psphotBlendFit.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+bool psphotBlendFitSetSource(pmSource *source, psImage *IDimage, float sigSigma);
 
 // for now, let's store the detections on the readout->analysis for each readout
@@ -92 +93 @@
     float skySig = psMetadataLookupF32(&status, recipe, "SKY_SIG");
     assert (status && isfinite(skySig) && skySig > 0);
+
+# if (0)
+    // **** test block : generate an ID image where pixels are set based the source models (flux > 0.1 peak && flux > 2.0 skySig)
+    psImage *IDimage = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_S32);
+    psImageInit (IDimage, 0);
+
+    // start with the currently known moments (Mxx, Myy, Mxy) and generate a window image
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+        psphotBlendFitSetSource (source, IDimage, skySig);
+    }
+    psphotSaveImage (NULL, IDimage, "idimage.fits");
+# endif
 
     // Define source fitting parameters for extended source fits
@@ -241 +255 @@
         pmSource *source = sources->data[i];
 
-# if (0)
-# define TEST_X 34
-# define TEST_Y 28
-   
+        int TEST_ON = false;
+# if (1)
+# define TEST_X 653
+# define TEST_Y 466
         if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
             fprintf (stderr, "test object\n");
-        }
-
+            //  psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
+            TEST_ON = true;
+        }
 # undef TEST_X
 # undef TEST_Y
@@ -258 +273 @@
         if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
         if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+
+        // skip saturated stars modeled with a radial profile 
+        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
 
         // skip DBL second sources (ie, added by psphotFitBlob
@@ -306 +324 @@
         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
             if (psphotFitBlob (readout, source, newSources, psf, fitOptions, maskVal, markVal)) {
+                if (TEST_ON) {
+                    psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
+                }
                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
                 Next ++;
@@ -312 +333 @@
         } else {
             if (psphotFitBlend (readout, source, psf, fitOptions, maskVal, markVal)) {
+                if (TEST_ON) {
+                    psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
+                }
                 source->type = PM_SOURCE_TYPE_STAR;
                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
@@ -349 +373 @@
 }
 
+bool psphotBlendFitSetSource(pmSource *source, psImage *IDimage, float skySigma) {
+
+    if (!source) return false;
+    if (!source->peak) return false; // XXX how can we have a peak-less source?
+    if (!source->moments) return false;
+    if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
+    if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
+    if (source->mode2 == PM_SOURCE_MODE2_MATCHED) return false;
+    psAssert(IDimage, "need a window");
+
+    if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) {
+        // find the radius at which we hit something like 0.1*SATURATION
+        return false;
+    }
+
+    if (!isfinite(source->moments->Mrf) || source->moments->Mrf < 0 ) return false;
+
+    // we have a source with moments Mx, My, Mxx, Mxy, Myy.  we just need to define a Gaussian that has 
+    // these values (and peak of 1.0)
+
+    int Nx = IDimage->numCols;
+    int Ny = IDimage->numRows;
+
+    float Xo = source->moments->Mx;
+    float Yo = source->moments->My;
+
+    psEllipseMoments moments;
+    moments.x2 = source->moments->Mxx;
+    moments.y2 = source->moments->Myy;
+    moments.xy = source->moments->Mxy;
+
+    psEllipseAxes axes = psEllipseMomentsToAxes(moments, 20.0);
+    if (! (isfinite(axes.major) && isfinite(axes.minor) && isfinite(axes.theta)) ) {
+        fprintf (stderr, "invalid axes found id: %4d major: %f minor: %f theta: %f\n", source->id, axes.major, axes.minor, axes.theta);
+        return false;
+    }
+
+    // Use 1st radial moment as size, not sigma.  Why this factor of 0.5 ?
+    float scale = 0.5 * source->moments->Mrf / axes.major;
+    axes.major *= scale;
+    axes.minor *= scale;
+
+    psEllipseShape shape = psEllipseAxesToShape(axes);
+    if (! (isfinite(shape.sx) && isfinite(shape.sy) && isfinite(shape.sxy)) ) {
+        fprintf (stderr, "invalid shape found id: %d sx: %f sy %f sxy: %f\n", source->id, shape.sx, shape.sy, shape.sxy);
+        return false;
+    }
+
+    float Smajor = axes.major;
+
+    // we want to go to 2.15 sigma (0.1 Io)
+    int minX = PS_MIN(PS_MAX(Xo - 2.15*Smajor, 0), Nx - 1);
+    int maxX = PS_MIN(PS_MAX(Xo + 2.15*Smajor, 0), Nx - 1);
+    int minY = PS_MIN(PS_MAX(Yo - 2.15*Smajor, 0), Ny - 1);
+    int maxY = PS_MIN(PS_MAX(Yo + 2.15*Smajor, 0), Ny - 1);
+
+    float rMxx = 0.5 / PS_SQR(shape.sx);
+    float rMyy = 0.5 / PS_SQR(shape.sy);
+    float Sxy = -1. * shape.sxy;    // factor of -1 is included to match the previous window function
+                                    // implementation. XXX: Is this correct?
+
+    int ID = source->id;
+    float Io = source->peak->rawFlux;
+    for (int iy = minY; iy < maxY; iy++) {
+        for (int ix = minX; ix < maxX; ix++) {
+
+            float dX = (ix + 0.5 - Xo);
+            float dY = (iy + 0.5 - Yo);
+
+            float z = rMxx * PS_SQR(dX) + rMyy * PS_SQR(dY) + Sxy*dX*dY;
+
+            float f = Io*exp(-z);
+            
+            if (f < 2.0*skySigma) continue;
+            IDimage->data.S32[iy][ix] = ID;
+        }
+    }
+
+    return true;
+}

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;
+}
+

trunk/psphot/src/psphotExtendedSourceAnalysis.c

@@ -218 +218 @@
         if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
 
+        // skip saturated stars modeled with a radial profile 
+        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
+
         // optionally allow non-extended objects to get petrosians as well
         if (!doPetroStars) {

trunk/psphot/src/psphotExtendedSourceFits.c

@@ -317 +317 @@
         if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
         if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+
+        // skip saturated stars modeled with a radial profile 
+        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
 
         // XXX this should use peak?

trunk/psphot/src/psphotFindDetections.c

@@ -120 +120 @@
     if (useFootprints) {
         if (replaceSourcesForFootprints) {
+            // subtract the noise for all sources including satstars
             psphotAddOrSubNoiseReadout(config, view, filerule, index, recipe, false);
             psphotReplaceAllSourcesReadout (config, view, filerule, index, recipe, false);
+
+            // add in the satstars
+            psphotAddOrSubSatstarsReadout (config, view, filerule, index, recipe, true);
+
             psFree (significance);
             significance = psphotSignificanceImage (readout, recipe, maskVal);
+
+            // display the significance image
+            psphotVisualShowSignificance (significance->variance, 0.98*threshold, 1.02*threshold);
         }
 
@@ -130 +138 @@
         if (replaceSourcesForFootprints) {
             psphotRemoveAllSourcesReadout (config, view, filerule, index, recipe, false);
+
+            // subtract the satstars
+            psphotAddOrSubSatstarsReadout (config, view, filerule, index, recipe, false);
         }
     }

trunk/psphot/src/psphotFitSourcesLinear.c

@@ -222 +222 @@
         if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
         if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+
+        // skip saturated stars modeled with a radial profile 
+        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
 
         // do not include CRs in the full ensemble fit

trunk/psphot/src/psphotGuessModels.c

@@ -169 +169 @@
         source->tmpFlags |= PM_SOURCE_TMPF_MODEL_GUESS;
 
-        // skip non-astronomical objects (very likely defects)
-        if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+        // skip non-astronomical objects (very likely defects) and satstars with profiles subtracted
+        if (source->mode  & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
         if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
         if (source->type == PM_SOURCE_TYPE_SATURATED) continue;

trunk/psphot/src/psphotKronIterate.c

@@ -331 +331 @@
             }
 
+            // skip saturated stars modeled with a radial profile 
+            if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
+
             // replace object in image
             bool reSubtract = false;

trunk/psphot/src/psphotMagnitudes.c

@@ -182 +182 @@
         if (saveTest) {
             psphotSaveImage(NULL, testImage, "test.image.1.fits");
+        }
+
+        // satstars modeled as a radial profile need special handling
+        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) {
+            psphotSatstarPhotometry (source);
+            continue;
         }
 

trunk/psphot/src/psphotRadialApertures.c

@@ -240 +240 @@
         if (source->mode & PM_SOURCE_MODE_DEFECT) continue;
 
+        // skip saturated stars modeled with a radial profile 
+        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
+
         // XXX measure radial apertures even for saturated stars
         // if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;

trunk/psphot/src/psphotRadialProfileWings.c

@@ -190 +190 @@
         if (!source->moments) continue;
 
+        // skip saturated stars modeled with a radial profile 
+        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
+
+        // skip non-detected sources matched from other images
+        if (source->mode2 & PM_SOURCE_MODE2_MATCHED) return true; // skip matched sources (no signal)
+
+        // XXX unclear if I should run this analysis on both 1st and 2nd pass for 1st pass objects,
+        // or just use the 1st pass value.  I think I should just use the 1st pass value, so skip
+        // any that have already been assigned
+        if (isfinite(source->skyRadius)) return true;
+
         // replace object in image
         bool reSubtract = false;
@@ -226 +237 @@
 
     // psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
-
-    // XXX unclear if I should run this analysis on both 1st and 2nd pass for 1st pass objects,
-    // or just use the 1st pass value.  I think I should just use the 1st pass value, so skip
-    // any that have already been assigned
-    if (isfinite(source->skyRadius)) return true;
-
-    if (source->mode2 & PM_SOURCE_MODE2_MATCHED) return true; // skip matched sources (no signal)
 
     // radii will be MIN_RADIUS to MAX_RADIUS in NN log steps:

trunk/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
@@ -149 +143 @@
     }
 
+    // find and subtract radial profile models for saturated stars (XXX change name eventually)
+    if (!psphotDeblendSatstars (config, view, filerule)) {
+        psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
+        return psphotReadoutCleanup (config, view, filerule);
+    }
+
     // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
     if (!psphotImageQuality (config, view, filerule)) { // pass 1
@@ -221 +221 @@
     // NOTE: this block performs the 2nd pass low-significance PSF detection stage
     { 
-        // add noise for subtracted objects
+        // add noise for subtracted objects & subtracted saturated stars
         psphotAddNoise (config, view, filerule); // pass 1 (detections->allSources)
 
@@ -253 +253 @@
         // merge the newly selected sources into the existing list
         // NOTE: merge OLD and NEW
-        // XXX check on free of sources...
         psphotMergeSources (config, view, filerule); // (detections->newSources + detections->allSources -> detections->allSources)
 
@@ -294 +293 @@
         // merge the newly selected sources into the existing list
         // NOTE: merge OLD and NEW
-        // XXX check on free of sources...
         psphotMergeSources (config, view, filerule); // (detections->newSources + detections->allSources -> detections->allSources)
 

trunk/psphot/src/psphotSourceSize.c

@@ -498 +498 @@
             continue;
         }
+
+        // skip saturated stars modeled with a radial profile 
+        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
 
         // we are classifying by moments and PSF_MAG - KRON_MAG

trunk/psphot/src/psphotSourceStats.c

@@ -447 +447 @@
         }
 
+        // skip saturated stars modeled with a radial profile (this probably never happens, since it is set after)
+        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
+
         if (!(source->peak->type == PM_PEAK_SUSPECT_SATURATION)) {
             // measure basic source moments (no S/N clipping on input pixels)
@@ -521 +524 @@
     float PSF_SN_LIM = 2.0*psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM");
 
-
+    // XXX this will cause an error in the vector length is > 8
     # define NSIGMA 8
     // moved to config file
@@ -538 +541 @@
             sigma[i] = sigmavec->data.F32[i];
         }
+        assert (sigmavec->n <= 8);
         nsigma = sigmavec->n;
     } else {

trunk/psphot/src/psphotStackMatchPSFs.c

@@ -46 +46 @@
     }
 
-    // loop over the available readouts (ignore chisq image)psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Source Stats ---");
+    // loop over the available readouts (ignore chisq image)
     for (int i = 0; i < num; i++) { 
         if (!psphotStackMatchPSFsReadout (config, view, options, i)) {

trunk/psphot/src/psphotVisual.c

@@ -262 +262 @@
 }
 
+static psImage *posImage = NULL;
+static psImage *delImage = NULL;
+
 bool psphotVisualShowImage (pmReadout *readout) {
 
@@ -277 +280 @@
     psphotVisualScaleImage (kapa, readout->variance, readout->mask, "variance", 1.0, 1);
     psphotVisualScaleImage (kapa, readout->image, readout->mask, "image", sqrt(factor), 0);
+
+    if (posImage == NULL) {
+        posImage = psImageCopy (NULL, readout->image, PS_TYPE_F32);
+    }
 
     pmVisualAskUser(NULL);
@@ -1139 +1146 @@
     // after displaying (as an image) the psf stars, we cycle throught them and display their
     // radial profiles:
-    psphotVisualPlotRadialProfiles (recipe, sources);
+    psphotVisualPlotRadialProfiles (recipe, sources, PM_SOURCE_MODE_PSFSTAR);
 
     return true;
@@ -1273 +1280 @@
 }
 
-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 +1292 @@
 }
 
-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);
@@ -1302 +1311 @@
     int ng = 0;
     int nb = 0;
-    float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS] - source->pixels->col0;
-    float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - source->pixels->row0;
+
+    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]) {
+            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) ;
@@ -1322 +1341 @@
     }
 
-    // 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(300, PS_TYPE_F32);
-    psVector *fmaj = psVectorAlloc(300, PS_TYPE_F32);
-    psVector *fmin = psVectorAlloc(300, 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);
-
     KapaInitGraph (&graphdata);
 
@@ -1392 +1348 @@
     // 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;
@@ -1418 +1374 @@
     KapaPlotVector (myKapa, nb, fb->data.F32, "y");
 
-    graphdata.color = KapaColorByName ("blue");
-    graphdata.ptype = 0;
-    graphdata.size = 0.0;
-    graphdata.style = 0;
-    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
-    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
-    KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
-    plotline (myKapa, &graphdata, 0.0, logHM, 30.0, logHM);
-    plotline (myKapa, &graphdata, 0.5*FWHM_MINOR, 0.0, 0.5*FWHM_MINOR, 5.0);
-    graphdata.ltype = 1;
-    plotline (myKapa, &graphdata, 0.5*MOMENTS_MINOR, 0.0, 0.5*MOMENTS_MINOR, 5.0);
-    graphdata.ltype = 0;
-        
-    graphdata.color = KapaColorByName ("green");
-    graphdata.ptype = 0;
-    graphdata.size = 0.0;
-    graphdata.style = 0;
-    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
-    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
-    KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
-    plotline (myKapa, &graphdata, 0.5*FWHM_MAJOR, 0.0, 0.5*FWHM_MAJOR, 5.0);
-    graphdata.ltype = 1;
-    plotline (myKapa, &graphdata, 0.5*MOMENTS_MAJOR, 0.0, 0.5*MOMENTS_MAJOR, 5.0);
-    graphdata.ltype = 0;
-        
-    for (int i = 0; i < rmod->n; i++) {
-        rmod->data.F32[i] = log10(rmod->data.F32[i]);
-    }
-
     // ** loglog **
     KapaSelectSection (myKapa, "loglog");
@@ -1452 +1379 @@
     // 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;
@@ -1479 +1406 @@
     KapaPlotVector (myKapa, nb, fb->data.F32, "y");
 
-    graphdata.color = KapaColorByName ("blue");
-    graphdata.ptype = 0;
-    graphdata.size = 0.0;
-    graphdata.style = 0;
-    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
-    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
-    KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
-
-    graphdata.color = KapaColorByName ("green");
-    graphdata.ptype = 0;
-    graphdata.size = 0.0;
-    graphdata.style = 0;
-    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
-    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
-    KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
-
-    psFree (rmod);
-    psFree (fmin);
-    psFree (fmaj);
+    if (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 (myKapa, "linlog");
+        KapaGetGraphData (myKapa, &graphdata);
+        graphdata.color = KapaColorByName ("blue");
+        graphdata.ptype = 0;
+        graphdata.size = 0.0;
+        graphdata.style = 0;
+        KapaPrepPlot   (myKapa, rmod->n, &graphdata);
+        KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
+        KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
+        plotline (myKapa, &graphdata, graphdata.xmin, logHM, graphdata.xmax, logHM);
+        plotline (myKapa, &graphdata, 0.5*FWHM_MINOR, graphdata.ymin, 0.5*FWHM_MINOR, graphdata.ymax);
+        graphdata.ltype = 1;
+        plotline (myKapa, &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   (myKapa, rmod->n, &graphdata);
+        KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
+        KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
+        plotline (myKapa, &graphdata, 0.5*FWHM_MAJOR, graphdata.ymin, 0.5*FWHM_MAJOR, graphdata.ymax);
+        graphdata.ltype = 1;
+        plotline (myKapa, &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 (myKapa, "loglog");
+        KapaGetGraphData (myKapa, &graphdata);
+        graphdata.color = KapaColorByName ("blue");
+        graphdata.ptype = 0;
+        graphdata.size = 0.0;
+        graphdata.style = 0;
+        KapaPrepPlot   (myKapa, rmod->n, &graphdata);
+        KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
+        KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
+
+        graphdata.color = KapaColorByName ("green");
+        graphdata.ptype = 0;
+        graphdata.size = 0.0;
+        graphdata.style = 0;
+        KapaPrepPlot   (myKapa, rmod->n, &graphdata);
+        KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
+        KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
+
+        psFree (rmod);
+        psFree (fmin);
+        psFree (fmaj);
+    }
 
     psFree (rg);
@@ -1508 +1535 @@
 }
 
-bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources) {
+bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources, pmSourceMode showmode) {
 
     KapaSection section;  // put the positive profile in one and the residuals in another?
@@ -1549 +1576 @@
 
         pmSource *source = sources->data[i];
-        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
-
-        psphotVisualPlotRadialProfile (myKapa, source, maskVal);
+        if (!(source->mode & showmode)) continue;
+
+        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:
@@ -2540 +2587 @@
     }
 
-    if (reshow) {
+    if (false && reshow) {
         psphotVisualShowMask (myKapa, readout->mask, "mask", 2);
         psphotVisualScaleImage (myKapa, readout->variance, readout->mask, "variance", 1.0, 1);
     }
-    psphotVisualScaleImage (myKapa, readout->image, readout->mask, "resid", sqrt(factor), 0);
+
+    if (posImage) {
+        delImage = (psImage *) psBinaryOp(delImage, posImage, "-", readout->image);
+        psphotVisualScaleImage (myKapa, posImage, readout->mask, "posimage", sqrt(factor), 0);
+        psphotVisualScaleImage (myKapa, delImage, readout->mask, "delimage", sqrt(factor), 2);
+    }
+
+    psphotVisualScaleImage (myKapa, readout->image, readout->mask, "resid", sqrt(factor), 1);
 
     pmVisualAskUser(NULL);