Changeset 34360

Timestamp:

Aug 26, 2012, 1:05:27 PM (14 years ago)

Author:

eugene

Message:

visual adjustment to examine saturated stars

Location:

branches/eam_branches/ipp-20120805/psphot/src

Files:

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

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

@@ -282 +282 @@
 bool            psphotVisualShowPSFModel (pmReadout *readout, pmPSF *psf);
 bool            psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal);
-bool            psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources);
+bool            psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources, pmSourceMode showmode);
 bool            psphotVisualShowFlags (psArray *sources);
 bool            psphotVisualShowSourceSize (pmReadout *readout, psArray *sources);

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

@@ -78 +78 @@
     psVector *index = psVectorSortIndex (NULL, SN);
     // this results in an index of increasing SN
+
+    psphotVisualPlotRadialProfiles (recipe, sources, PM_SOURCE_MODE_SATSTAR);
 
     // examine sources in decreasing SN order
@@ -219 +221 @@
     return true;
 }
+
+bool psphotDeblendSatstarsReadoutOld (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex) {
+
+    int N;
+    pmSource *source;
+    bool status;
+
+    psTimerStart ("psphot.deblend.sat");
+
+    int Nblend = 0;
+    float SAT_MIN_RADIUS = 5.0;
+
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, fileIndex); // File of interest
+    psAssert (file, "missing file?");
+
+    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+    psAssert (readout, "missing readout?");
+
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detections, "missing detections?");
+
+    psArray *sources = detections->newSources;
+    psAssert (sources, "missing sources?");
+
+    if (!sources->n) {
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping satstar blend");
+        return true;
+    }
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    pmCell *cell = readout->parent;
+
+    float SATURATION = NAN;
+
+    // do not completely trust the values in the header...
+    float CELL_SATURATION = psMetadataLookupF32 (&status, cell->concepts, "CELL.SATURATION");
+    float MIN_SATURATION = psMetadataLookupF32 (&status, recipe, "DEBLEND_MIN_SATURATION");
+    if (!status || !isfinite(MIN_SATURATION)) {
+        MIN_SATURATION = 40000.0;
+    }
+    if (!isfinite(CELL_SATURATION)) {
+        SATURATION = MIN_SATURATION;
+    } else {
+        SATURATION = PS_MAX(MIN_SATURATION, CELL_SATURATION);
+    }
+    float SAT_TEST_LEVEL = 0.5*SATURATION;
+
+    // we need sources spatially-sorted to find overlaps
+    sources = psArraySort (sources, pmSourceSortByY);
+
+    // source analysis is done in peak order (brightest first)
+    // we use an index for this so the spatial sorting is kept
+    psVector *SN = psVectorAlloc (sources->n, PS_DATA_F32);
+    for (int i = 0; i < SN->n; i++) {
+        source = sources->data[i];
+        SN->data.F32[i] = source->peak->rawFlux;
+    }
+    psVector *index = psVectorSortIndex (NULL, SN);
+    // this results in an index of increasing SN
+
+    // examine sources in decreasing SN order
+    for (int i = sources->n - 1; i >= 0; i--) {
+        N = index->data.U32[i];
+        source = sources->data[N];
+
+        // XXX filter? if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
+        if (source->mode & PM_SOURCE_MODE_BLEND) continue;
+        if (source->peak->rawFlux < SAT_TEST_LEVEL) continue;
+
+        // save these for reference below
+        int xPeak = source->peak->x;
+        int yPeak = source->peak->y;
+
+        // generate a basic contour (set of x,y coordinates at-or-below flux level)
+        psArray *contour = pmSourceContour (source->pixels, xPeak, yPeak, SAT_TEST_LEVEL);
+        if (contour == NULL) continue;
+
+        // contour consists of a set of X,Y coords giving the boundary
+        psVector *xVec = contour->data[0];
+        psVector *yVec = contour->data[1];
+        if (xVec->n < 5) {
+            psFree(contour);
+            continue;
+        }
+
+        // find the center of the contour (let's just use mid[x,y])
+        int xMin = xVec->data.F32[0];
+        int xMax = xVec->data.F32[0];
+        int yMin = yVec->data.F32[0];
+        int yMax = yVec->data.F32[0];
+        for (int j = 0; j < xVec->n; j++) {
+            xMin = PS_MIN (xMin, xVec->data.F32[j]);
+            xMax = PS_MAX (xMax, xVec->data.F32[j]);
+            yMin = PS_MIN (yMin, yVec->data.F32[j]);
+            yMax = PS_MAX (yMax, yVec->data.F32[j]);
+        }       
+        int xCenter = 0.5*(xMin + xMax);
+        int yCenter = 0.5*(yMin + yMax);
+        psFree (contour);
+
+        psAssert (xCenter >= source->pixels->col0, "invalid shift in object center");
+        psAssert (xCenter <  source->pixels->col0 + source->pixels->numCols, "invalid shift in object center");
+        psAssert (yCenter >= source->pixels->row0, "invalid shift in object center");
+        psAssert (yCenter <  source->pixels->row0 + source->pixels->numRows, "invalid shift in object center");
+
+        // reset the peak for this source to the value of the center pixel
+        source->peak->x = xCenter;
+        source->peak->y = yCenter;
+        source->peak->xf = xCenter;
+        source->peak->yf = yCenter;
+        
+        // temporary array for overlapping objects we find
+        psArray *overlap = psArrayAllocEmpty (100);
+
+        // search backwards for overlapping sources
+        for (int j = N - 1; j >= 0; j--) {
+            pmSource *testSource = sources->data[j];
+            if (testSource->peak->x <  source->pixels->col0) continue;
+            if (testSource->peak->x >= source->pixels->col0 + source->pixels->numCols) continue;
+            if (testSource->peak->y <  source->pixels->row0) break;
+            if (testSource->peak->y >= source->pixels->row0 + source->pixels->numRows) {
+                fprintf (stderr, "warning: invalid condition\n");
+                continue;
+            }
+            psArrayAdd (overlap, 100, testSource);
+        }
+
+        // search forwards for overlapping sources
+        for (int j = N + 1; j < sources->n; j++) {
+            pmSource *testSource = sources->data[j];
+            if (testSource->peak->x <  source->pixels->col0) continue;
+            if (testSource->peak->x >= source->pixels->col0 + source->pixels->numCols) continue;
+            if (testSource->peak->y <  source->pixels->row0) {
+                fprintf (stderr, "warning: invalid condition\n");
+                continue;
+            }
+            if (testSource->peak->y >= source->pixels->row0 + source->pixels->numRows) break;
+            psArrayAdd (overlap, 100, testSource);
+        }
+
+        if (overlap->n == 0) {
+            psFree (overlap);
+            continue;
+        }
+
+        // now find the contour which is at 0.5*SAT_TEST_LEVEL (xPeak, yPeak is valid high pixel)
+        contour = pmSourceContour (source->pixels, xPeak, yPeak, 0.5*SAT_TEST_LEVEL);
+        if (contour == NULL) {
+            psFree (overlap);
+            continue; 
+        }
+        
+        // any peaks within this contour should be dropped (mark as blend, drop after loop is done)
+        // also drop any peaks which are too close to this peak
+        psVector *xv = contour->data[0];
+        psVector *yv = contour->data[1];
+        for (int k = 0; k < overlap->n; k++) {
+            pmSource *testSource = overlap->data[k];
+            float radius = hypot((testSource->peak->x - xCenter), (testSource->peak->y - yCenter));
+            if (radius < SAT_MIN_RADIUS) {
+                testSource->mode |= PM_SOURCE_MODE_BLEND;
+                Nblend ++;
+                continue;
+            }
+            for (int j = 0; j < xv->n; j+=2) {
+                if (fabs(yv->data.F32[j] - testSource->peak->y) > 0.5) continue;
+                if (xv->data.F32[j+0] > testSource->peak->x) break;
+                if (xv->data.F32[j+1] < testSource->peak->x) break;
+                testSource->mode |= PM_SOURCE_MODE_BLEND;
+                Nblend ++;
+                j = xv->n; // skip rest of contour
+            }
+        }
+        psFree (overlap);
+        psFree (contour);
+    }
+
+    // drop the sources marked as BLEND
+    for (int i = 0; i < sources->n;) {
+        pmSource *source = sources->data[i];
+
+        if (!(source->mode & PM_SOURCE_MODE_BLEND)) {
+            i++;
+            continue;
+        }
+
+        // shuffle the remaining sources forward
+        for (int j = i; j < sources->n - 1; j++) {
+            sources->data[j] = sources->data[j+1];
+        }
+        psFree (source);
+        sources->n --;
+    }
+
+    psFree (SN);
+    psFree (index);
+
+    psLogMsg ("psphot", PS_LOG_INFO, "found %d satstar blend peaks, leaving %ld sources: %f sec\n", Nblend, sources->n, psTimerMark ("psphot.deblend.sat"));
+    return true;
+}

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

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

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

@@ -1139 +1139 @@
     // 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;
@@ -1302 +1302 @@
     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;
+        Yo = source->moments->My;
+    }
+
     for (int iy = 0; iy < source->pixels->numRows; iy++) {
         for (int ix = 0; ix < source->pixels->numCols; ix++) {
@@ -1322 +1332 @@
     }
 
-    // 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);
 
@@ -1417 +1364 @@
     KapaPlotVector (myKapa, nb, rb->data.F32, "x");
     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 **
@@ -1479 +1397 @@
     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(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);
+
+        // ** 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 +1526 @@
 }
 
-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 +1567 @@
 
         pmSource *source = sources->data[i];
-        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
+        if (!(source->mode & showmode)) continue;
 
         psphotVisualPlotRadialProfile (myKapa, source, maskVal);