Index: /branches/eam_branch_20080926/psphot/src/psphotVisual.c
===================================================================
--- /branches/eam_branch_20080926/psphot/src/psphotVisual.c	(revision 19861)
+++ /branches/eam_branch_20080926/psphot/src/psphotVisual.c	(revision 19862)
@@ -1,3 +1,15 @@
 # include "psphotInternal.h"
+bool psphotVisualPlotMoments (pmConfig *config, const pmFPAview *view, psArray *sources);
+
+// this function displays representative images as the psphot analysis progresses:
+// 0 : image, 1 : weight
+// 0 : backsub, 1 : weight, 2 : backgnd
+// 0 : backsub, 1 : weight, 2 : signif
+// (overlay peaks on images)
+// (overlay footprints on images)
+// (overlay moments on images)
+// (overlay rough class on images) 
+// 0 : backsub, 1 : psfpos, 2: psfsub
+// 0 : backsub, 1 : lin_resid, 2: psfsub
 
 # if (HAVE_KAPA)
@@ -9,69 +21,1072 @@
 static bool isVisual = false;
 static int kapa = -1;
+static int kapa2 = -1;
+static int kapa3 = -1;
 
 bool psphotSetVisual (bool mode) {
 
-  isVisual = mode;
-  return true;
-}
-
-bool psphotVisualShowImage (pmReadout *readout) {
-
-  KiiImage image;
-  KapaImageData data;
-  Coords coords;
-
-  if (!isVisual) return true;
-
-  if (kapa == -1) {
-    kapa = pmKapaOpen (true, "psphot");
+    isVisual = mode;
+    return true;
+}
+
+bool psphotVisualScaleImage (int kapaFD, psImage *inImage, const char *name, int channel) {
+
+    KiiImage image;
+    KapaImageData data;
+    Coords coords;
+
+    strcpy (coords.ctype, "RA---TAN");
+
+    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
+    if (!psImageBackground(stats, NULL, inImage, NULL, 0, rng)) {
+	fprintf (stderr, "failed to get background values\n");
+	return false;
+    }
+
+    image.data2d = inImage->data.F32;
+    image.Nx = inImage->numCols;
+    image.Ny = inImage->numRows;
+
+    strcpy (data.name, name);
+    strcpy (data.file, name);
+    data.zero = stats->robustMedian - stats->robustStdev;
+    data.range = 5*stats->robustStdev;
+    data.logflux = 0;
+  
+    KiiSetChannel (kapaFD, channel);
+    KiiNewPicture2D (kapaFD, &image, &data, &coords);
+
+    psFree (stats);
+    psFree (rng);
+
+    return true;
+}
+
+bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max) {
+
+    KiiImage image;
+    KapaImageData data;
+    Coords coords;
+
+    strcpy (coords.ctype, "RA---TAN");
+
+    image.data2d = inImage->data.F32;
+    image.Nx = inImage->numCols;
+    image.Ny = inImage->numRows;
+
+    strcpy (data.name, name);
+    strcpy (data.file, name);
+    data.zero = min;
+    data.range = max - min;
+    data.logflux = 0;
+  
+    KiiSetChannel (kapaFD, channel);
+    KiiNewPicture2D (kapaFD, &image, &data, &coords);
+
+    return true;
+}
+
+bool psphotVisualShowImage (pmConfig *config, pmReadout *readout) {
+
+    if (!isVisual) return true;
+
     if (kapa == -1) {
-      fprintf (stderr, "failure to open kapa; visual mode disabled\n");
-      isVisual = false;
-      return false;
-    }
-  }  
-
-  // XXX get image stats
-  image.data2d = readout->image->data.F32;
-  image.Nx = readout->image->numCols;
-  image.Ny = readout->image->numRows;
-
-  strcpy (data.name, "image");
-  strcpy (data.file, "image");
-  data.zero = 0.0;
-  data.range = 1000.0;
-  data.logflux = 0;
+        kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
+	if (kapa == -1) {
+	    fprintf (stderr, "failure to open kapa; visual mode disabled\n");
+	    isVisual = false;
+	    return false;
+	}
+    }  
+
+    // psphotVisualShowMask (kapa, readout->mask, "mask", 2);
+    psphotVisualScaleImage (kapa, readout->weight, "weight", 1);
+    psphotVisualScaleImage (kapa, readout->image, "image", 0);
+
+    // pause and wait for user input:
+    // continue, save (provide name), ??
+    char key[10];
+    fprintf (stdout, "[c]ontinue? ");
+    if (!fgets(key, 8, stdin)) {
+	psWarning("Unable to read option");
+    }
+    return true;
+}
+
+bool psphotVisualShowBackground (pmConfig *config, const pmFPAview *view, pmReadout *readout) {
+
+    pmReadout *backgnd;
+
+    if (!isVisual) return true;
+
+    if (kapa == -1) {
+        kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
+	if (kapa == -1) {
+	    fprintf (stderr, "failure to open kapa; visual mode disabled\n");
+	    isVisual = false;
+	    return false;
+	}
+    }  
+
+    bool status = false;
+    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKGND");
+
+    if (file->mode == PM_FPA_MODE_INTERNAL) {
+	backgnd = file->readout;
+    } else {
+	backgnd = pmFPAviewThisReadout (view, file->fpa);
+    }
+
+    psphotVisualScaleImage (kapa, backgnd->image, "backgnd", 2);
+    psphotVisualScaleImage (kapa, readout->image, "backsub", 0);
+
+    // pause and wait for user input:
+    // continue, save (provide name), ??
+    char key[10];
+    fprintf (stdout, "[c]ontinue? ");
+    if (!fgets(key, 8, stdin)) {
+	psWarning("Unable to read option");
+    }
+    return true;
+}
+
+bool psphotVisualShowSignificance (psImage *image) {
+
+    if (!isVisual) return true;
+
+    if (kapa == -1) {
+        kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
+	if (kapa == -1) {
+	    fprintf (stderr, "failure to open kapa; visual mode disabled\n");
+	    isVisual = false;
+	    return false;
+	}
+    }  
+
+    // XXX test: image->data.F32[10][10] = 10000;
+    psphotVisualRangeImage (kapa, image, "signif", 2, -1.0, 25.0*25.0);
+
+    // pause and wait for user input:
+    // continue, save (provide name), ??
+    char key[10];
+    fprintf (stdout, "[c]ontinue? ");
+    if (!fgets(key, 8, stdin)) {
+	psWarning("Unable to read option");
+    }
+    return true;
+}
+
+bool psphotVisualShowPeaks (pmConfig *config, const pmFPAview *view, pmDetections *detections) {
+
+    int Noverlay;
+    KiiOverlay *overlay;
   
-  KiiSetChannel (kapa, 0);
-  KiiNewPicture2D (kapa, &image, &data, &coords);
+    if (!isVisual) return true;
+
+    if (kapa == -1) {
+	fprintf (stderr, "kapa not opened, skipping\n");
+	return false;
+    }
+
+    psArray *peaks = detections->peaks;
+
+    // note: this uses the Ohana allocation tools:
+    // ALLOCATE (overlay, KiiOverlay, 3*peaks->n + 1);
+    ALLOCATE (overlay, KiiOverlay, peaks->n);
+
+    Noverlay = 0;
+    for (int i = 0; i < peaks->n; i++) {
+
+	pmPeak *peak = peaks->data[i];
+	if (peak == NULL) continue;
+
+	overlay[Noverlay].type = KII_OVERLAY_BOX;
+	overlay[Noverlay].x = peak->xf;
+	overlay[Noverlay].y = peak->yf;
+	overlay[Noverlay].dx = 2.0;
+	overlay[Noverlay].dy = 2.0;
+	overlay[Noverlay].angle = 0.0;
+	overlay[Noverlay].text = NULL;
+	Noverlay ++;
+
+# if (0)
+	overlay[Noverlay].type = KII_OVERLAY_BOX;
+	overlay[Noverlay].x = peak->x;
+	overlay[Noverlay].y = peak->y;
+	overlay[Noverlay].dx = 1.0;
+	overlay[Noverlay].dy = 1.0;
+	overlay[Noverlay].angle = 0.0;
+	overlay[Noverlay].text = NULL;
+	Noverlay ++;
+
+	overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
+	overlay[Noverlay].x = peak->xf;
+	overlay[Noverlay].y = peak->yf;
+	overlay[Noverlay].dx = 2.0;
+	overlay[Noverlay].dy = 2.0;
+	overlay[Noverlay].angle = 0.0;
+	overlay[Noverlay].text = NULL;
+	Noverlay ++;
+# endif
+    }
+
+# if (0)
+    overlay[Noverlay].type = KII_OVERLAY_BOX;
+    overlay[Noverlay].x = 10.0;
+    overlay[Noverlay].y = 10.0;
+    overlay[Noverlay].dx = 0.5;
+    overlay[Noverlay].dy = 0.5;
+    overlay[Noverlay].angle = 0.0;
+    overlay[Noverlay].text = NULL;
+    Noverlay ++;
+# endif
+
+    KiiLoadOverlay (kapa, overlay, Noverlay, "red");
+    FREE (overlay);
+
+    // pause and wait for user input:
+    // continue, save (provide name), ??
+    char key[10];
+    fprintf (stdout, "[c]ontinue? ");
+    if (!fgets(key, 8, stdin)) {
+	psWarning("Unable to read option");
+    }
+    return true;
+}
+
+// XXX this is just wrong: outlines entire image
+bool psphotVisualShowFootprints (pmConfig *config, const pmFPAview *view, pmDetections *detections) {
+
+    int Noverlay;
+    KiiOverlay *overlay;
   
-  // display the weight
-  image.data2d = readout->weight->data.F32;
-  image.Nx = readout->weight->numCols;
-  image.Ny = readout->weight->numRows;
-
-  strcpy (data.name, "weight");
-  strcpy (data.file, "weight");
-  data.zero = 0.0;
-  data.range = 1000.0;
-  data.logflux = 0;
+    if (!isVisual) return true;
+
+    if (kapa == -1) {
+	fprintf (stderr, "kapa not opened, skipping\n");
+	return false;
+    }
+
+    psArray *footprints = detections->footprints;
+    if (!footprints) return true;
+
+    // note: this uses the Ohana allocation tools:
+    int NOVERLAY = footprints->n;
+    ALLOCATE (overlay, KiiOverlay, NOVERLAY);
+
+    Noverlay = 0;
+    for (int i = 0; i < footprints->n; i++) {
+
+	pmSpan *span = NULL;
+
+	pmFootprint *footprint = footprints->data[i];
+	if (footprint == NULL) continue;
+	if (footprint->spans == NULL) continue;
+	if (footprint->spans->n < 1) continue;
+
+	// draw the top
+	span = footprint->spans->data[0];
+	overlay[Noverlay].type = KII_OVERLAY_LINE;
+	overlay[Noverlay].x = span->x0;
+	overlay[Noverlay].y = span->y;
+	overlay[Noverlay].dx = span->x1 - span->x0;
+	overlay[Noverlay].dy = 0;
+	overlay[Noverlay].angle = 0.0;
+	overlay[Noverlay].text = NULL;
+	Noverlay ++;
+	CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
+
+	int ys = span->y;
+	int x0s = span->x0;
+	int x1s = span->x1;
+
+	// draw the outer span edges
+	for (int j = 1; j < footprint->spans->n; j++) {
+	    pmSpan *span1 = footprint->spans->data[j];
+
+	    int ye = span1->y;
+	    int x0e = span1->x0;
+	    int x1e = span1->x1;
+
+	    // we cannot have two discontinuous spans on the top or bottom, right? (no, probably not right)
+	    // find all of the spans in this row and generate x0e, x01:
+	    for (int k = j + 1; k < footprint->spans->n; k++) {
+		pmSpan *span2 = footprint->spans->data[k];
+		if (span2->y > span1->y) break;
+		x0e = PS_MIN (x0e, span2->x0);
+		x1e = PS_MAX (x1e, span2->x1);
+		j++;
+	    }
+
+	    overlay[Noverlay].type = KII_OVERLAY_LINE;
+	    overlay[Noverlay].x = x0s;
+	    overlay[Noverlay].y = ys;
+	    overlay[Noverlay].dx = x0e - x0s;
+	    overlay[Noverlay].dy = ye - ys;
+	    overlay[Noverlay].angle = 0.0;
+	    overlay[Noverlay].text = NULL;
+	    Noverlay ++;
+	    CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
+
+	    overlay[Noverlay].type = KII_OVERLAY_LINE;
+	    overlay[Noverlay].x = x1s;
+	    overlay[Noverlay].y = ys;
+	    overlay[Noverlay].dx = x1e - x1s;
+	    overlay[Noverlay].dy = ye - ys;
+	    overlay[Noverlay].angle = 0.0;
+	    overlay[Noverlay].text = NULL;
+	    Noverlay ++;
+	    CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
+
+	    ys = ye;
+	    x0s = x0e;
+	    x1s = x1e;
+	}
+
+	// draw the bottom
+	span = footprint->spans->data[footprint->spans->n - 1];
+	overlay[Noverlay].type = KII_OVERLAY_LINE;
+	overlay[Noverlay].x = span->x0;
+	overlay[Noverlay].y = span->y;
+	overlay[Noverlay].dx = span->x1 - span->x0;
+	overlay[Noverlay].dy = 0;
+	overlay[Noverlay].angle = 0.0;
+	overlay[Noverlay].text = NULL;
+	Noverlay ++;
+	CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
+    }
+
+    KiiLoadOverlay (kapa, overlay, Noverlay, "blue");
+    FREE (overlay);
+
+    // pause and wait for user input:
+    // continue, save (provide name), ??
+    char key[10];
+    fprintf (stdout, "[c]ontinue? ");
+    if (!fgets(key, 8, stdin)) {
+	psWarning("Unable to read option");
+    }
+    return true;
+}
+
+bool psphotVisualShowMoments (pmConfig *config, const pmFPAview *view, psArray *sources) {
+
+    int Noverlay;
+    KiiOverlay *overlay;
   
-  KiiSetChannel (kapa, 0);
-  KiiNewPicture2D (kapa, &image, &data, &coords);
+    psEllipseMoments emoments;
+    psEllipseAxes axes;
+
+    if (!isVisual) return true;
+
+    if (kapa == -1) {
+	fprintf (stderr, "kapa not opened, skipping\n");
+	return false;
+    }
+
+    // note: this uses the Ohana allocation tools:
+    ALLOCATE (overlay, KiiOverlay, sources->n);
+
+    Noverlay = 0;
+    for (int i = 0; i < sources->n; i++) {
+
+	pmSource *source = sources->data[i];
+	if (source == NULL) continue;
+
+	pmMoments *moments = source->moments;
+	if (moments == NULL) continue;
+
+	overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
+	overlay[Noverlay].x = moments->x;
+	overlay[Noverlay].y = moments->y;
+
+	emoments.x2 = moments->Sx;
+	emoments.y2 = moments->Sy;
+	emoments.xy = moments->Sxy;
+
+	axes = psEllipseMomentsToAxes (emoments, 20.0);
+
+	overlay[Noverlay].dx = 2.0*axes.major;
+	overlay[Noverlay].dy = 2.0*axes.minor;
+	overlay[Noverlay].angle = axes.theta * PS_DEG_RAD;
+	overlay[Noverlay].text = NULL;
+	Noverlay ++;
+    }
+
+    KiiLoadOverlay (kapa, overlay, Noverlay, "yellow");
+    FREE (overlay);
+
+    // pause and wait for user input:
+    // continue, save (provide name), ??
+    char key[10];
+    fprintf (stdout, "[c]ontinue? ");
+    if (!fgets(key, 8, stdin)) {
+	psWarning("Unable to read option");
+    }
+
+    return true;
+}
+
+bool psphotVisualPlotMoments (pmConfig *config, const pmFPAview *view, psArray *sources) {
+
+    bool status;
+    Graphdata graphdata;
+
+    if (!isVisual) return true;
+
+    if (kapa3 == -1) {
+        kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots");
+	if (kapa3 == -1) {
+	    fprintf (stderr, "failure to open kapa; visual mode disabled\n");
+	    isVisual = false;
+	    return false;
+	}
+    }  
+
+    // select the current recipe
+    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
+    if (!recipe) {
+	fprintf (stderr, "missing recipe, skipping moments plot\n");
+        return false;
+    }
+
+    KapaClearPlots (kapa3);
+    KapaInitGraph (&graphdata);
+
+    float psfX  = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.X");
+    float psfY  = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.Y");
+    float psfdX = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DX");
+    float psfdY = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DY");
+    float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
+    float Rx = psfdX * PSF_CLUMP_NSIGMA;
+    float Ry = psfdY * PSF_CLUMP_NSIGMA;
+
+    // examine sources to set data range
+    graphdata.xmin = -0.05;
+    graphdata.ymin = -0.05;
+    graphdata.xmax = 2.0*psfX;
+    graphdata.ymax = 2.0*psfY;
+    KapaSetLimits (kapa3, &graphdata);
+
+    KapaSetFont (kapa3, "helvetica", 14);
+    KapaBox (kapa3, &graphdata);
+    KapaSendLabel (kapa3, "&ss&h_x| (pixels)", KAPA_LABEL_XM);
+    KapaSendLabel (kapa3, "&ss&h_y| (pixels)", KAPA_LABEL_YM);
+
+    psVector *xBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *yBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *xFaint  = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *yFaint  = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+
+    // construct the vectors
+    int nB = 0;
+    int nF = 0;
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+        if (source->moments == NULL)
+            continue;
+
+        xFaint->data.F32[nF] = source->moments->Sx;
+        yFaint->data.F32[nF] = source->moments->Sy;
+        nF++;
+
+        // XXX make this a user-defined cutoff
+        if (source->moments->SN < 50)
+            continue;
+
+        xBright->data.F32[nB] = source->moments->Sx;
+        yBright->data.F32[nB] = source->moments->Sy;
+        nB++;
+    }
+    xFaint->n = nF;
+    yFaint->n = nF;
+
+    xBright->n = nB;
+    yBright->n = nB;
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 0;
+    graphdata.size = 0.3;
+    graphdata.style = 2;
+    KapaPrepPlot (kapa3, nF, &graphdata);
+    KapaPlotVector (kapa3, nF, xFaint->data.F32, "x");
+    KapaPlotVector (kapa3, nF, yFaint->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("red");
+    graphdata.ptype = 2;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+    KapaPrepPlot (kapa3, nB, &graphdata);
+    KapaPlotVector (kapa3, nB, xBright->data.F32, "x");
+    KapaPlotVector (kapa3, nB, yBright->data.F32, "y");
+
+    // draw a circle centered on psfX,Y with size of the psf limit
+    psVector *xLimit  = psVectorAlloc (120, PS_TYPE_F32);
+    psVector *yLimit  = psVectorAlloc (120, PS_TYPE_F32);
+    for (int i = 0; i < xLimit->n; i++) {
+	xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX;
+	yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY;
+    }
+    graphdata.color = KapaColorByName ("blue");
+    graphdata.style = 0;
+    KapaPrepPlot (kapa3, xLimit->n, &graphdata);
+    KapaPlotVector (kapa3, xLimit->n, xLimit->data.F32, "x");
+    KapaPlotVector (kapa3, yLimit->n, yLimit->data.F32, "y");
+    psFree (xLimit);
+    psFree (yLimit);
+
+# if (0)
+    // *** make a histogram of the source counts in the x and y directions 
+    psHistogram *nX = psHistogramAlloc (graphdata.xmin, graphdata.xmax, 50.0);
+    psHistogram *nY = psHistogramAlloc (graphdata.ymin, graphdata.ymax, 50.0);
+    psVectorHistogram (nX, xFaint, NULL, NULL, 0);
+    psVectorHistogram (nY, yFaint, NULL, NULL, 0);
+    psVector *dX = psVectorAlloc (nX->nums->n, PS_TYPE_F32);
+    psVector *vX = psVectorAlloc (nX->nums->n, PS_TYPE_F32);
+    psVector *dY = psVectorAlloc (nY->nums->n, PS_TYPE_F32);
+    psVector *vY = psVectorAlloc (nY->nums->n, PS_TYPE_F32);
+    for (int i = 0; i < nX->nums->n; i++) {
+	dX->data.F32[i] = nX->nums->data.S32[i];
+	vX->data.F32[i] = 0.5*(nX->bounds->data.F32[i] + nX->bounds->data.F32[i+1]);
+    }
+    for (int i = 0; i < nY->nums->n; i++) {
+	dY->data.F32[i] = nY->nums->data.S32[i];
+	vY->data.F32[i] = 0.5*(nY->bounds->data.F32[i] + nY->bounds->data.F32[i+1]);
+    }
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 0;
+    graphdata.size = 0.0;
+    graphdata.style = 0;
+    KapaPrepPlot (kapa3, dX->n, &graphdata);
+    KapaPlotVector (kapa3, dX->n, dX->data.F32, "x");
+    KapaPlotVector (kapa3, vX->n, vX->data.F32, "y");
+
+    psFree (nX);
+    psFree (dX);
+    psFree (vX);
+
+    psFree (nY);
+    psFree (dY);
+    psFree (vY);
+# endif
+
+    psFree (xBright);
+    psFree (yBright);
+    psFree (xFaint);
+    psFree (yFaint);
+
+    // pause and wait for user input:
+    // continue, save (provide name), ??
+    char key[10];
+    fprintf (stdout, "[c]ontinue? ");
+    if (!fgets(key, 8, stdin)) {
+	psWarning("Unable to read option");
+    }
+    return true;
+}
+
+// assumes 'kapa' value is checked and set
+bool psphotVisualShowRoughClass_Single (pmConfig *config, const pmFPAview *view, psArray *sources, pmSourceType type, pmSourceMode mode, char *color) {
+
+    int Noverlay;
+    KiiOverlay *overlay;
   
-  // display the weight
-  image.data2d = readout->weight->data.F32;
-  image.Nx = readout->weight->numCols;
-  image.Ny = readout->weight->numRows;
-
-  strcpy (data.name, "mask");
-  strcpy (data.file, "maak");
-  data.zero = 0.0;
-  data.range = 1000.0;
-  data.logflux = 0;
-  
-  KiiSetChannel (kapa, 0);
-  KiiNewPicture2D (kapa, &image, &data, &coords);
-}
-
+    psEllipseMoments emoments;
+    psEllipseAxes axes;
+
+    // note: this uses the Ohana allocation tools:
+    ALLOCATE (overlay, KiiOverlay, sources->n);
+
+    Noverlay = 0;
+    for (int i = 0; i < sources->n; i++) {
+
+	pmSource *source = sources->data[i];
+	if (source == NULL) continue;
+
+	if (source->type != type) continue;
+	if (mode && !(source->mode & mode)) continue;
+
+	pmMoments *moments = source->moments;
+	if (moments == NULL) continue;
+
+	overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
+	overlay[Noverlay].x = moments->x;
+	overlay[Noverlay].y = moments->y;
+
+	emoments.x2 = moments->Sx;
+	emoments.y2 = moments->Sy;
+	emoments.xy = moments->Sxy;
+
+	axes = psEllipseMomentsToAxes (emoments, 20.0);
+
+	overlay[Noverlay].dx = 2.0*axes.major;
+	overlay[Noverlay].dy = 2.0*axes.minor;
+	overlay[Noverlay].angle = axes.theta * PS_DEG_RAD;
+	overlay[Noverlay].text = NULL;
+	Noverlay ++;
+    }
+
+    KiiLoadOverlay (kapa, overlay, Noverlay, color);
+    FREE (overlay);
+
+    return true;
+}
+
+bool psphotVisualShowRoughClass (pmConfig *config, const pmFPAview *view, psArray *sources) {
+
+    psphotVisualPlotMoments (config, view, sources);
+
+    if (!isVisual) return true;
+
+    if (kapa == -1) {
+	fprintf (stderr, "kapa not opened, skipping\n");
+	return false;
+    }
+
+    KiiEraseOverlay (kapa, "yellow"); // moments
+
+    psphotVisualShowRoughClass_Single (config, view, sources, PM_SOURCE_TYPE_STAR, 0, "red");
+    psphotVisualShowRoughClass_Single (config, view, sources, PM_SOURCE_TYPE_EXTENDED, 0, "blue");
+    psphotVisualShowRoughClass_Single (config, view, sources, PM_SOURCE_TYPE_DEFECT, 0, "blue");
+    psphotVisualShowRoughClass_Single (config, view, sources, PM_SOURCE_TYPE_SATURATED, 0, "red");
+    psphotVisualShowRoughClass_Single (config, view, sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_PSFSTAR, "yellow");
+    psphotVisualShowRoughClass_Single (config, view, sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_SATSTAR, "green");
+
+    // pause and wait for user input:
+    // continue, save (provide name), ??
+    char key[10];
+    fprintf (stdout, "red: STAR or SAT AREA; blue: EXTENDED or DEFECT; green: SATSTAR; yellow: PSFSTAR\n");
+    fprintf (stdout, "[c]ontinue? ");
+    if (!fgets(key, 8, stdin)) {
+	psWarning("Unable to read option");
+    }
+    return true;
+}
+
+bool psphotVisualShowPSFStars (pmConfig *config, const pmFPAview *view, pmPSF *psf, psArray *sources) {
+
+    bool status;
+
+    if (!isVisual) return true;
+
+    if (kapa2 == -1) {
+        kapa2 = KapaOpenNamedSocket ("kapa", "psphot:psfstars");
+	if (kapa2 == -1) {
+	    fprintf (stderr, "failure to open kapa; visual mode disabled\n");
+	    isVisual = false;
+	    return false;
+	}
+    }  
+
+    // select the current recipe
+    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
+    if (!recipe) {
+        psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE);
+        return false;
+    }
+
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
+    assert (maskVal);
+
+    // the source images are written to an image 10x the size of a PSF object
+    // float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
+    // PS_ASSERT (status, false);
+
+    int DX = 21;
+    int DY = 21;
+
+    // examine PSF sources in S/N order (brightest first)
+    sources = psArraySort (sources, pmSourceSortBySN);
+
+    // counters to track the size of the image and area used in a row
+    int dX = 0;				// starting corner of next box
+    int dY = 0;				// height of row so far
+    int NX = 20*DX;			// full width of output image
+    int NY = 0;				// total height of output image
+
+    // first, examine the PSF stars:
+    // - determine bounding boxes for summary image
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *source = sources->data[i];
+
+	bool keep = false;
+        keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR);
+	if (!keep) continue;
+
+	// how does this subimage get placed into the output image?
+	// DX = source->pixels->numCols
+	// DY = source->pixels->numRows
+
+	if (dX + DX > NX) {
+	    // too wide for the rest of this row
+	    if (dX == 0) {
+		// alone on this row
+		NY += DY;
+		dX = 0;
+		dY = 0;
+	    } else {
+		// start the next row
+		NY += dY;
+		dX = DX;
+		dY = DY;
+	    }
+	} else {
+	    // extend this row
+	    dX += DX;
+	    dY = PS_MAX (dY, DY);
+	}
+    }
+    NY += DY;
+
+    // allocate output image
+    psImage *outpos = psImageAlloc (NX, NY, PS_TYPE_F32);
+    psImage *outsub = psImageAlloc (NX, NY, PS_TYPE_F32);
+
+    int Xo = 0;				// starting corner of next box
+    int Yo = 0;				// starting corner of next box
+    dY = 0;				// height of row so far
+
+    int nPSF = 0;
+
+    // next, examine the PSF stars:
+    // - create output image array
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *source = sources->data[i];
+
+	bool keep = false;
+        if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
+	    nPSF ++;
+	    keep = true;
+	}
+	if (!keep) continue;
+
+	if (Xo + DX > NX) {
+	    // too wide for the rest of this row
+	    if (Xo == 0) {
+		// place source alone on this row
+		psphotAddWithTest (source, true, maskVal); // replace source if subtracted
+		psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
+
+		psphotSubWithTest (source, false, maskVal); // remove source (force)
+		psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
+		psphotSetState (source, false, maskVal); // reset source Add/Sub state to recorded
+
+		Yo += DY;
+		Xo = 0;
+		dY = 0;
+	    } else {
+		// start the next row
+		Yo += dY;
+		Xo = 0;
+
+		psphotAddWithTest (source, true, maskVal); // replace source if subtracted
+		psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
+
+		psphotSubWithTest (source, false, maskVal); // remove source (force)
+		psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
+		psphotSetState (source, false, maskVal); // replace source (has been subtracted)
+
+		Xo = DX;
+		dY = DY;
+	    }
+	} else {
+	    // extend this row
+	    psphotAddWithTest (source, true, maskVal); // replace source if subtracted
+	    psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
+
+	    psphotSubWithTest (source, false, maskVal); // remove source (force)
+	    psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
+	    psphotSetState (source, false, maskVal); // replace source (has been subtracted)
+
+	    Xo += DX;
+	    dY = PS_MAX (dY, DY);
+	}
+    }
+
+    psphotVisualScaleImage (kapa2, outpos, "psfpos", 0);
+    psphotVisualScaleImage (kapa2, outsub, "psfsub", 1);
+
+    // pause and wait for user input:
+    // continue, save (provide name), ??
+    char key[10];
+    fprintf (stdout, "[c]ontinue? ");
+    if (!fgets(key, 8, stdin)) {
+	psWarning("Unable to read option");
+    }
+
+    psFree (outpos);
+    psFree (outsub);
+    return true;
+}
+
+bool psphotVisualShowSatStars (pmConfig *config, const pmFPAview *view, pmPSF *psf, psArray *sources) {
+
+    bool status;
+
+    if (!isVisual) return true;
+
+    if (kapa2 == -1) {
+        kapa2 = KapaOpenNamedSocket ("kapa", "psphot:images");
+	if (kapa2 == -1) {
+	    fprintf (stderr, "failure to open kapa; visual mode disabled\n");
+	    isVisual = false;
+	    return false;
+	}
+    }  
+
+    // select the current recipe
+    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
+    if (!recipe) {
+        psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE);
+        return false;
+    }
+
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
+    assert (maskVal);
+
+    // the source images are written to an image 10x the size of a PSF object
+    // float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
+    // PS_ASSERT (status, false);
+
+    int DX = 41;
+    int DY = 41;
+
+    // examine PSF sources in S/N order (brightest first)
+    sources = psArraySort (sources, pmSourceSortBySN);
+
+    // counters to track the size of the image and area used in a row
+    int dX = 0;				// starting corner of next box
+    int dY = 0;				// height of row so far
+    int NX = 10*DX;			// full width of output image
+    int NY = 0;				// total height of output image
+
+    // first, examine the PSF and SAT stars:
+    // - determine bounding boxes for summary image
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *source = sources->data[i];
+
+	bool keep = false;
+        keep |= (source->mode & PM_SOURCE_MODE_SATSTAR);
+	if (!keep) continue;
+
+	// how does this subimage get placed into the output image?
+	// DX = source->pixels->numCols
+	// DY = source->pixels->numRows
+
+	if (dX + DX > NX) {
+	    // too wide for the rest of this row
+	    if (dX == 0) {
+		// alone on this row
+		NY += DY;
+		dX = 0;
+		dY = 0;
+	    } else {
+		// start the next row
+		NY += dY;
+		dX = DX;
+		dY = DY;
+	    }
+	} else {
+	    // extend this row
+	    dX += DX;
+	    dY = PS_MAX (dY, DY);
+	}
+    }
+    NY += DY;
+
+    // allocate output image
+    psImage *outsat = psImageAlloc (NX, NY, PS_TYPE_F32);
+
+    int Xo = 0;				// starting corner of next box
+    int Yo = 0;				// starting corner of next box
+    dY = 0;				// height of row so far
+
+    int nSAT = 0;
+
+    // next, examine the SAT stars:
+    // - create output image array
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *source = sources->data[i];
+
+	bool keep = false;
+        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
+	    nSAT ++;
+	    keep = true;
+	}	    
+	if (!keep) continue;
+
+	if (Xo + DX > NX) {
+	    // too wide for the rest of this row
+	    if (Xo == 0) {
+		// place source alone on this row
+		psphotAddWithTest (source, true, maskVal); // replace source if subtracted
+		psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
+		psphotSetState (source, true, maskVal); // reset source Add/Sub state to recorded
+
+		Yo += DY;
+		Xo = 0;
+		dY = 0;
+	    } else {
+		// start the next row
+		Yo += dY;
+		Xo = 0;
+		psphotAddWithTest (source, true, maskVal); // replace source if subtracted
+		psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
+		psphotSetState (source, true, maskVal); // replace source (has been subtracted)
+
+		Xo = DX;
+		dY = DY;
+	    }
+	} else {
+	    // extend this row
+	    psphotAddWithTest (source, true, maskVal); // replace source if subtracted
+	    psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
+	    psphotSetState (source, true, maskVal); // replace source (has been subtracted)
+
+	    Xo += DX;
+	    dY = PS_MAX (dY, DY);
+	}
+    }
+
+    psphotVisualScaleImage (kapa2, outsat, "satstar", 2);
+
+    // pause and wait for user input:
+    // continue, save (provide name), ??
+    char key[10];
+    fprintf (stdout, "[c]ontinue? ");
+    if (!fgets(key, 8, stdin)) {
+	psWarning("Unable to read option");
+    }
+
+    psFree (outsat);
+    return true;
+}
+
+bool psphotVisualShowPSFModel (pmConfig *config, pmReadout *readout, pmPSF *psf) {
+
+    if (!isVisual) return true;
+
+    if (kapa2 == -1) {
+        kapa2 = KapaOpenNamedSocket ("kapa", "psphot:psfstars");
+	if (kapa2 == -1) {
+	    fprintf (stderr, "failure to open kapa; visual mode disabled\n");
+	    isVisual = false;
+	    return false;
+	}
+    }  
+
+    int DX = 64;
+    int DY = 64;
+
+    psImage *psfMosaic = psImageAlloc (3*DX, 3*DY, PS_TYPE_F32);
+    psImageInit (psfMosaic, 0.0);
+    
+    pmModel *modelRef = pmModelAlloc(psf->type);
+
+# if (1)
+    // generate a fake model at each of the 3x3 image grid positions
+    for (int x = -1; x <= +1; x ++) {
+	for (int y = -1; y <= +1; y ++) {
+	    // use the center of the center pixel of the image
+	    float xc = (0.5 + 0.45*x)*readout->image->numCols + readout->image->col0;
+	    float yc = (0.5 + 0.45*y)*readout->image->numRows + readout->image->row0;
+
+	    // assign the x and y coords to the image center
+	    // create an object with center intensity of 1000
+	    modelRef->params->data.F32[PM_PAR_SKY] = 0;
+	    modelRef->params->data.F32[PM_PAR_I0] = 1000;
+	    modelRef->params->data.F32[PM_PAR_XPOS] = xc;
+	    modelRef->params->data.F32[PM_PAR_YPOS] = yc;
+    
+	    // create modelPSF from this model
+	    pmModel *model = pmModelFromPSF (modelRef, psf);
+
+	    // place the reference object in the image center
+	    // no need to mask the source here
+	    // XXX should we measure this for the analytical model only or the full model?
+	    pmModelAddWithOffset (psfMosaic, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_CENTER, 0, -x*DX, -y*DY);
+	}
+    }
+# else
+    { 
+	// use the center of the center pixel of the image
+	float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5;
+	float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5;
+
+	// assign the x and y coords to the image center
+	// create an object with center intensity of 1000
+	modelRef->params->data.F32[PM_PAR_SKY] = 0;
+	modelRef->params->data.F32[PM_PAR_I0] = 1000;
+	modelRef->params->data.F32[PM_PAR_XPOS] = xc;
+	modelRef->params->data.F32[PM_PAR_YPOS] = yc;
+    
+	// create modelPSF from this model
+	pmModel *model = pmModelFromPSF (modelRef, psf);
+
+	// place the reference object in the image center
+	// no need to mask the source here
+	// XXX should we measure this for the analytical model only or the full model?
+	pmModelAddWithOffset (psfMosaic, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_CENTER, 0, 32.0, 16.0);
+    }
+# endif
+
+    psImage *psfLogFlux = (psImage *) psUnaryOp (NULL, psfMosaic, "log");
+    psphotVisualRangeImage (kapa2, psfLogFlux, "psf_mosaic", 1, -2.0, 3.0);
+
+    // pause and wait for user input:
+    // continue, save (provide name), ??
+    char key[10];
+    fprintf (stdout, "[c]ontinue? ");
+    if (!fgets(key, 8, stdin)) {
+	psWarning("Unable to read option");
+    }
+    return true;
+}
+
+bool psphotVisualShowLinearFit (pmConfig *config, pmReadout *readout) {
+
+    if (!isVisual) return true;
+
+    if (kapa == -1) {
+        kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
+	if (kapa == -1) {
+	    fprintf (stderr, "failure to open kapa; visual mode disabled\n");
+	    isVisual = false;
+	    return false;
+	}
+    }  
+
+    psphotVisualScaleImage (kapa, readout->image, "lin_resid", 1);
+
+    // pause and wait for user input:
+    // continue, save (provide name), ??
+    char key[10];
+    fprintf (stdout, "[c]ontinue? ");
+    if (!fgets(key, 8, stdin)) {
+	psWarning("Unable to read option");
+    }
+    return true;
+}
+
+# else
+
+bool psphotSetVisual (bool mode){}
+bool psphotVisualShowImage (pmConfig *config, pmReadout *readout){}
+bool psphotVisualShowBackground (pmConfig *config, const pmFPAview *view, pmReadout *readout){}
+bool psphotVisualShowSignificance (psImage *image){}
+bool psphotVisualShowPeaks (pmConfig *config, const pmFPAview *view, pmDetections *detections){}
+bool psphotVisualShowFootprints (pmConfig *config, const pmFPAview *view, pmDetections *detections){}
+bool psphotVisualShowMoments (pmConfig *config, const pmFPAview *view, psArray *sources){}
+bool psphotVisualShowRoughClass (pmConfig *config, const pmFPAview *view, psArray *sources){}
+bool psphotVisualShowPSF (pmConfig *config, const pmFPAview *view, pmPSF *psf, psArray *sources){}
+bool psphotVisualShowLinearFit (pmConfig *config, pmReadout *readout){}
+
+# endif