# include "psphotInternal.h"

// 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)
# include <kapa.h>

// functions used to visualize the analysis as it goes
// these are invoked by the -visual options

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 psphotVisualShowMask (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.Nx = inImage->numCols;
    image.Ny = inImage->numRows;

    ALLOCATE (image.data2d, float *, image.Ny);
    for (int iy = 0; iy < image.Ny; iy++) {
	ALLOCATE (image.data2d[iy], float, image.Nx);
	for (int ix = 0; ix < image.Nx; ix++) {
	    image.data2d[iy][ix] = inImage->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix];
	}
    }

    strcpy (data.name, name);
    strcpy (data.file, name);
    data.zero = -1;
    data.range = 32;
    data.logflux = 0;
  
    KiiSetChannel (kapaFD, channel);
    KiiNewPicture2D (kapaFD, &image, &data, &coords);

    for (int iy = 0; iy < image.Ny; iy++) {
	free (image.data2d[iy]);
    }
    free (image.data2d);

    psFree (stats);
    psFree (rng);

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

    // 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 (pmDetections *detections) {

    int Noverlay;
    KiiOverlay *overlay;
  
    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;
}

bool psphotVisualShowFootprints (pmDetections *detections) {

    int Noverlay;
    KiiOverlay *overlay;
  
    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 (psArray *sources) {

    int Noverlay;
    KiiOverlay *overlay;
  
    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->Mx;
	overlay[Noverlay].y = moments->My;

	emoments.x2 = moments->Mxx;
	emoments.xy = moments->Mxy;
	emoments.y2 = moments->Myy;

	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;  // XXXXXXXX the axes angle is negative to display of object on kapa
	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 (psMetadata *recipe, 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;
	}
    }  

    KapaClearPlots (kapa3);
    KapaInitGraph (&graphdata);

    // there are N regions: use the first (guaranteed to exist) to get the overall limits
    psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, "PSF.CLUMP.REGION.000");
    
    float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
    float psfY  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");

    // 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->Mxx;
        yFaint->data.F32[nF] = source->moments->Myy;
        nF++;

        // XXX make this a user-defined cutoff
        if (source->moments->SN < 50)
            continue;

        xBright->data.F32[nB] = source->moments->Mxx;
        yBright->data.F32[nB] = source->moments->Myy;
        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");

    // XXX draw N circles to outline the clumps
    { 
	// 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);

	int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
	float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");

	graphdata.color = KapaColorByName ("blue");
	graphdata.style = 0;

	for (int n = 0; n < nRegions; n++) {

	    char regionName[64];
	    snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
	    psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
    
	    float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
	    float psfY  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
	    float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
	    float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
	    float Rx = psfdX * PSF_CLUMP_NSIGMA;
	    float Ry = psfdY * PSF_CLUMP_NSIGMA;

	    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;
	    }
	    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 (psArray *sources, pmSourceType type, pmSourceMode mode, char *color) {

    int Noverlay;
    KiiOverlay *overlay;
  
    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->Mx;
	overlay[Noverlay].y = moments->My;

	emoments.x2 = moments->Mxx;
	emoments.y2 = moments->Myy;
	emoments.xy = moments->Mxy;

	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 (psArray *sources) {

    if (!isVisual) return true;

    if (kapa == -1) {
	fprintf (stderr, "kapa not opened, skipping\n");
	return false;
    }

    KiiEraseOverlay (kapa, "yellow"); // moments

    psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_STAR, 0, "red");
    psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_EXTENDED, 0, "blue");
    psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_DEFECT, 0, "blue");
    psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_SATURATED, 0, "red");
    psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_PSFSTAR, "yellow");
    psphotVisualShowRoughClass_Single (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 psphotVisualShowPSFModel (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 (5*DX, 5*DY, PS_TYPE_F32);
    psImageInit (psfMosaic, 0.0);
    
    psImage *funMosaic = psImageAlloc (5*DX, 5*DY, PS_TYPE_F32);
    psImageInit (funMosaic, 0.0);

    psImage *resMosaic = psImageAlloc (5*DX, 5*DY, PS_TYPE_F32);
    psImageInit (resMosaic, 0.0);

    pmModel *modelRef = pmModelAlloc(psf->type);

    // generate a fake model at each of the 3x3 image grid positions
    for (int x = -2; x <= +2; x ++) {
	for (int y = -2; y <= +2; y ++) {
	    // use the center of the center pixel of the image
	    float xc = (int)((0.5 + 0.225*x)*readout->image->numCols) + readout->image->col0 + 0.5;
	    float yc = (int)((0.5 + 0.225*y)*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);
	    if (!model) continue;

	    // 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);
	    pmModelAddWithOffset (funMosaic, NULL, model, PM_MODEL_OP_FUNC | PM_MODEL_OP_CENTER, 0, x*DX, y*DY);
	    pmModelAddWithOffset (resMosaic, NULL, model, PM_MODEL_OP_RES0 | PM_MODEL_OP_RES1 | PM_MODEL_OP_CENTER, 0, x*DX, y*DY);
	    psFree (model);
	}
    }

    psImage *psfLogFlux = (psImage *) psUnaryOp (NULL, psfMosaic, "log");
    psphotVisualRangeImage (kapa2, psfLogFlux, "psf_mosaic",    0, -2.0, 3.0);
    psphotVisualRangeImage (kapa2, funMosaic, "psf_analytical", 1, -10.0, 100.0);
    psphotVisualRangeImage (kapa2, resMosaic, "psf_residual",   2, -10.0, 100.0);

    psFree (psfMosaic);
    psFree (funMosaic);
    psFree (resMosaic);
    psFree (psfLogFlux);
    psFree (modelRef);

    // 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 psphotVisualShowPSFStars (psMetadata *recipe, 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;
	}
    }  

    // 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);

    // 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);
    psImageInit (outpos, 0.0);
    psImageInit (outsub, 0.0);

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

    psphotVisualRangeImage (kapa2, outpos, "psfpos", 0, -0.05, 0.95);
    psphotVisualRangeImage (kapa2, outsub, "psfsub", 1, -0.05, 0.95);

    // 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);


    // after displaying (as an image) the psf stars, we cycle throught them and display their
    // radial profiles:
    psphotVisualPlotRadialProfiles (recipe, sources);

    return true;
}

bool psphotVisualShowSatStars (psMetadata *recipe, 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;
	}
    }  

    // 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);

    // 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);
    psImageInit (outsat, 0.0);

    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 psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal) {

    Graphdata graphdata;

    bool state = !(source->mode & PM_SOURCE_MODE_SUBTRACTED);
    psphotAddWithTest (source, true, maskVal); // replace source if subtracted

    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 *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 = source->modelPSF->params->data.F32[PM_PAR_XPOS] - source->pixels->col0;
    float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - 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]) {
		// 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] = log10(source->pixels->data.F32[iy][ix]);
		ng++;
	    }
	}
    }
  
    // reset source Add/Sub state to recorded
    psphotSetState (source, state, maskVal); 

    KapaInitGraph (&graphdata);

    // ** linlog **
    KapaSelectSection (myKapa, "linlog");

    // examine sources to set data range
    graphdata.xmin =  -0.05;
    graphdata.xmax = +30.05;
    graphdata.ymin = -0.05;
    graphdata.ymax = +5.05;
    KapaSetLimits (myKapa, &graphdata);
  
    KapaSetFont (myKapa, "helvetica", 14);
    KapaBox (myKapa, &graphdata);
    KapaSendLabel (myKapa, "radius (pixels)", KAPA_LABEL_XM);
    KapaSendLabel (myKapa, "log flux (counts)", KAPA_LABEL_YM);
	       
    graphdata.color = KapaColorByName ("black");
    graphdata.ptype = 2;
    graphdata.size = 0.5;
    graphdata.style = 2;
    KapaPrepPlot (myKapa, ng, &graphdata);
    KapaPlotVector (myKapa, ng, rg->data.F32, "x");
    KapaPlotVector (myKapa, ng, fg->data.F32, "y");

    graphdata.color = KapaColorByName ("red");
    graphdata.ptype = 0;
    graphdata.size = 0.3;
    graphdata.style = 2;
    KapaPrepPlot (myKapa, nb, &graphdata);
    KapaPlotVector (myKapa, nb, rb->data.F32, "x");
    KapaPlotVector (myKapa, nb, fb->data.F32, "y");
  
    // ** loglog **
    KapaSelectSection (myKapa, "loglog");

    // examine sources to set data range
    graphdata.xmin = -1.51;
    graphdata.xmax = +1.51;
    graphdata.ymin = -0.05;
    graphdata.ymax = +5.05;
    KapaSetLimits (myKapa, &graphdata);
  
    KapaSetFont (myKapa, "helvetica", 14);
    KapaBox (myKapa, &graphdata);
    KapaSendLabel (myKapa, "log radius (pixels)", KAPA_LABEL_XM);
    KapaSendLabel (myKapa, "log flux (counts)", KAPA_LABEL_YM);
	       
    graphdata.color = KapaColorByName ("black");
    graphdata.ptype = 2;
    graphdata.size = 0.5;
    graphdata.style = 2;
    KapaPrepPlot (myKapa, ng, &graphdata);
    KapaPlotVector (myKapa, ng, Rg->data.F32, "x");
    KapaPlotVector (myKapa, ng, fg->data.F32, "y");

    graphdata.color = KapaColorByName ("red");
    graphdata.ptype = 0;
    graphdata.size = 0.3;
    graphdata.style = 2;
    KapaPrepPlot (myKapa, nb, &graphdata);
    KapaPlotVector (myKapa, nb, Rb->data.F32, "x");
    KapaPlotVector (myKapa, nb, fb->data.F32, "y");

    psFree (rg);
    psFree (Rg);
    psFree (fg);
    psFree (rb);
    psFree (Rb);
    psFree (fb);
    return true;
}

bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources) {

    KapaSection section;  // put the positive profile in one and the residuals in another? 

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

    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    bool status;
    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    assert (maskVal);

    KapaClearPlots (kapa3);
    // first section : mag vs CR nSigma
    section.dx = 1.0;
    section.dy = 0.5;
    section.x = 0.0;
    section.y = 0.0;
    section.name = NULL;
    psStringAppend (&section.name, "linlog");
    KapaSetSection (kapa3, &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.name = NULL;
    psStringAppend (&section.name, "loglog");
    KapaSetSection (kapa3, &section);
    psFree (section.name);

    // loop over the PSF stars
    for (int i = 0; i < sources->n; i++) {

        pmSource *source = sources->data[i];
        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;

	psphotVisualPlotRadialProfile (kapa3, source, maskVal);

	// pause and wait for user input:
	// continue, save (provide name), ??
	char key[10];
	fprintf (stdout, "[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 (kapa3);
	}
	if (key[0] == 's') {
	    break;
	}
    }	    

    return true;
}

bool psphotVisualShowFlags (psArray *sources) {

    int NoverlayE, NOVERLAYE;
    int NoverlayO, NOVERLAYO;
    KiiOverlay *overlayE, *overlayO;
  
    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:
    NoverlayE = 0;
    NOVERLAYE = 100;
    ALLOCATE (overlayE, KiiOverlay, NOVERLAYE);

    NoverlayO = 0;
    NOVERLAYO = 100;
    ALLOCATE (overlayO, KiiOverlay, NOVERLAYO);

    for (int i = 0; i < sources->n; i++) {

	float Xo, Yo, Rmaj, Rmin, cs, sn;

	pmSource *source = sources->data[i];
	if (source == NULL) continue;

	pmMoments *moments = source->moments;
	if (0) { 
	    emoments.x2 = moments->Mxx;
	    emoments.y2 = moments->Myy;
	    emoments.xy = moments->Mxy;
	    Xo = moments->Mx;
	    Yo = moments->My;

	    axes = psEllipseMomentsToAxes (emoments, 20.0);
	    Rmaj = 2.0*axes.major;
	    Rmin = 2.0*axes.minor;
	    cs = cos(axes.theta);
	    sn = sin(axes.theta);
	} else {
	    Rmaj = Rmin = 5.0;
	    cs = 1.0;
	    sn = 0.0;
	    Xo = source->peak->xf;
	    Yo = source->peak->yf;
	}

	unsigned short int flagMask = 0x01;
	for (int j = 0; j < 8; j++) {
	    if (source->mode & flagMask) {
		overlayE[NoverlayE].type = KII_OVERLAY_LINE;
		overlayE[NoverlayE].x = Xo;
		overlayE[NoverlayE].y = Yo;

		float phi = j*M_PI/4.0;
		overlayE[NoverlayE].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn;
		overlayE[NoverlayE].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs;
		overlayE[NoverlayE].angle = 0;
		overlayE[NoverlayE].text = NULL;
		NoverlayE ++;
		CHECK_REALLOCATE (overlayE, KiiOverlay, NOVERLAYE, NoverlayE, 100);
	    }
	    flagMask <<= 1;

	    if (source->mode & flagMask) {
		overlayO[NoverlayO].type = KII_OVERLAY_LINE;
		overlayO[NoverlayO].x = Xo + 1;
		overlayO[NoverlayO].y = Yo;

		float phi = j*M_PI/4.0;
		overlayO[NoverlayO].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn;
		overlayO[NoverlayO].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs;
		overlayO[NoverlayO].angle = 0;
		overlayO[NoverlayO].text = NULL;
		NoverlayO ++;
		CHECK_REALLOCATE (overlayO, KiiOverlay, NOVERLAYO, NoverlayO, 100);
	    }
	    flagMask <<= 1;
	}
    }

    KiiLoadOverlay (kapa, overlayE, NoverlayE, "red");
    KiiLoadOverlay (kapa, overlayO, NoverlayO, "yellow");
    FREE (overlayE);
    FREE (overlayO);

    // pause and wait for user input:
    // continue, save (provide name), ??
    char key[10];
    fprintf (stdout, "even bits (0x0001, 0x0004, ... : red\n");
    fprintf (stdout, "odd bits (0x0002, 0x0008, ... : yellow\n");
    fprintf (stdout, "[c]ontinue? ");
    if (!fgets(key, 8, stdin)) {
	psWarning("Unable to read option");
    }

    return true;
}

bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources) {

    int Noverlay, NOVERLAY;
    KiiOverlay *overlay;

    if (!isVisual) return true;

    if (kapa == -1) {
	fprintf (stderr, "kapa not opened, skipping\n");
	return false;
    }

    // note: this uses the Ohana allocation tools:
    Noverlay = 0;
    NOVERLAY = 100;
    ALLOCATE (overlay, KiiOverlay, sources->n);

    // mark CRs with red boxes
    for (int i = 0; i < sources->n; i++) {

	pmSource *source = sources->data[i];
	if (source == NULL) continue;

	if (!(source->mode & PM_SOURCE_MODE_CR_LIMIT)) continue;

	overlay[Noverlay].type = KII_OVERLAY_BOX;
	overlay[Noverlay].x = source->peak->xf;
	overlay[Noverlay].y = source->peak->yf;

	overlay[Noverlay].dx = 4;
	overlay[Noverlay].dy = 4;
	overlay[Noverlay].angle = 0;
	overlay[Noverlay].text = NULL;
	Noverlay ++;
	CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    }
    KiiLoadOverlay (kapa, overlay, Noverlay, "red");


    Noverlay = 0;
    for (int i = 0; i < sources->n; i++) {

	pmSource *source = sources->data[i];
	if (source == NULL) continue;

	// mark EXTs with yellow circles
	if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;

	overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
	overlay[Noverlay].x = source->peak->xf;
	overlay[Noverlay].y = source->peak->yf;

	overlay[Noverlay].dx = 10;
	overlay[Noverlay].dy = 10;
	overlay[Noverlay].angle = 0;
	overlay[Noverlay].text = NULL;
	Noverlay ++;
	CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    }

    KiiLoadOverlay (kapa, overlay, Noverlay, "blue");
    FREE (overlay);

    psphotVisualShowMask (kapa, readout->mask, "mask", 2);

    // pause and wait for user input:
    // continue, save (provide name), ??
    char key[10];
    fprintf (stdout, "CR: 4pix red BOX; EXT: 10pix blue circle\n");
    fprintf (stdout, "[c]ontinue? ");
    if (!fgets(key, 8, stdin)) {
	psWarning("Unable to read option");
    }

    return true;
}

bool psphotVisualPlotSourceSize (psArray *sources) {

    Graphdata graphdata;
    KapaSection section;

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

    KapaClearPlots (kapa3);
    KapaInitGraph (&graphdata);

    // first section : mag vs CR nSigma
    section.dx = 1.0;
    section.dy = 0.5;
    section.x = 0.0;
    section.y = 0.0;
    section.name = NULL;
    psStringAppend (&section.name, "a1");
    KapaSetSection (kapa3, &section);
    psFree (section.name);

    psVector *x = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    psVector *y = psVectorAllocEmpty (sources->n, PS_TYPE_F32);

    graphdata.xmin = +32.0;
    graphdata.xmax = -32.0;
    graphdata.ymin = +32.0;
    graphdata.ymax = -32.0;

    // construct the plot vectors
    int n = 0;
    for (int i = 0; i < sources->n; i++) {
        pmSource *source = sources->data[i];
	if (!source) continue;
        if (source->type != PM_SOURCE_TYPE_STAR) continue;
	if (!isfinite (source->crNsigma)) continue;

        x->data.F32[n] = -2.5*log10(source->peak->flux);
        y->data.F32[n] = source->crNsigma;
        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
        graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
        graphdata.ymin = -0.5;
        graphdata.ymax = 10.0;

        n++;
    }
    x->n = y->n = n;

    float range;
    range = graphdata.xmax - graphdata.xmin;
    graphdata.xmax += 0.05*range;
    graphdata.xmin -= 0.05*range;

    // XXX set the plot range to match the image
    KapaSetLimits (kapa3, &graphdata);

    KapaSetFont (kapa3, "helvetica", 14);
    KapaBox (kapa3, &graphdata);
    KapaSendLabel (kapa3, "Peak as Mag", KAPA_LABEL_XM);
    KapaSendLabel (kapa3, "CR N Sigma", KAPA_LABEL_YM);

    graphdata.color = KapaColorByName ("black");
    graphdata.ptype = 2;
    graphdata.size = 0.5;
    graphdata.style = 2;
    KapaPrepPlot (kapa3, n, &graphdata);
    KapaPlotVector (kapa3, n, x->data.F32, "x");
    KapaPlotVector (kapa3, n, y->data.F32, "y");

    // second section : mag vs EXT nSigma
    section.dx = 1.0;
    section.dy = 0.5;
    section.x = 0.0;
    section.y = 0.5;
    section.name = NULL;
    psStringAppend (&section.name, "a2");
    KapaSetSection (kapa3, &section);
    psFree (section.name);

    graphdata.xmin = +32.0;
    graphdata.xmax = -32.0;
    graphdata.ymin = +32.0;
    graphdata.ymax = -32.0;

    // construct the plot vectors
    n = 0;
    for (int i = 0; i < sources->n; i++) {
        pmSource *source = sources->data[i];
	if (!source) continue;
        if (source->type != PM_SOURCE_TYPE_STAR) continue;
	if (!isfinite (source->extNsigma)) continue;

        x->data.F32[n] = -2.5*log10(source->peak->flux);
        y->data.F32[n] = source->extNsigma;
        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
        graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
        graphdata.ymin = -0.5;
        graphdata.ymax = 10.0;

        n++;
    }
    x->n = y->n = n;

    range = graphdata.xmax - graphdata.xmin;
    graphdata.xmax += 0.05*range;
    graphdata.xmin -= 0.05*range;

    // XXX set the plot range to match the image
    KapaSetLimits (kapa3, &graphdata);

    KapaSetFont (kapa3, "helvetica", 14);
    KapaBox (kapa3, &graphdata);
    KapaSendLabel (kapa3, "EXT N Sigma", KAPA_LABEL_YM);

    graphdata.color = KapaColorByName ("black");
    graphdata.ptype = 2;
    graphdata.size = 0.5;
    graphdata.style = 2;
    KapaPrepPlot (kapa3, n, &graphdata);
    KapaPlotVector (kapa3, n, x->data.F32, "x");
    KapaPlotVector (kapa3, n, y->data.F32, "y");

    psFree (x);
    psFree (y);

    // 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 psphotVisualShowResidualImage (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, "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;
}

bool psphotVisualPlotApResid (psArray *sources) {

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

    KapaClearPlots (kapa3);
    KapaInitGraph (&graphdata);

    psVector *x = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    psVector *y = psVectorAllocEmpty (sources->n, PS_TYPE_F32);

    graphdata.xmin = +32.0;
    graphdata.xmax = -32.0;
    graphdata.ymin = +32.0;
    graphdata.ymax = -32.0;

    // construct the plot vectors
    int n = 0;
    for (int i = 0; i < sources->n; i++) {
        pmSource *source = sources->data[i];
	if (!source) continue;
        if (source->type != PM_SOURCE_TYPE_STAR) continue;
	if (!isfinite (source->apMag)) continue;
	if (!isfinite (source->psfMag)) continue;

        x->data.F32[n] = source->psfMag;
        y->data.F32[n] = source->apMag - source->psfMag;
        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
        graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
        graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]);
        graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]);

        n++;
    }
    x->n = y->n = n;

    float range;
    range = graphdata.xmax - graphdata.xmin;
    graphdata.xmax += 0.05*range;
    graphdata.xmin -= 0.05*range;
    range = graphdata.ymax - graphdata.ymin;
    graphdata.ymax += 0.05*range;
    graphdata.ymin -= 0.05*range;

    // XXX set the plot range to match the image
    KapaSetLimits (kapa3, &graphdata);

    KapaSetFont (kapa3, "helvetica", 14);
    KapaBox (kapa3, &graphdata);
    KapaSendLabel (kapa3, "PSF Mag", KAPA_LABEL_XM);
    KapaSendLabel (kapa3, "Ap Mag - PSF Mag", KAPA_LABEL_YM);

    graphdata.color = KapaColorByName ("black");
    graphdata.ptype = 2;
    graphdata.size = 0.5;
    graphdata.style = 2;
    KapaPrepPlot (kapa3, n, &graphdata);
    KapaPlotVector (kapa3, n, x->data.F32, "x");
    KapaPlotVector (kapa3, n, y->data.F32, "y");

    psFree (x);
    psFree (y);

    // 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) { return true; }
bool psphotVisualShowBackground (pmConfig *config, const pmFPAview *view, pmReadout *readout) { return true; }
bool psphotVisualShowSignificance (psImage *image) { return true; }
bool psphotVisualShowPeaks (pmConfig *config, const pmFPAview *view, pmDetections *detections) { return true; }
bool psphotVisualShowFootprints (pmConfig *config, const pmFPAview *view, pmDetections *detections) { return true; }
bool psphotVisualShowMoments (pmConfig *config, const pmFPAview *view, psArray *sources) { return true; }
bool psphotVisualPlotMoments (pmConfig *config, const pmFPAview *view, psArray *sources) { return true; }
bool psphotVisualShowRoughClass (pmConfig *config, const pmFPAview *view, psArray *sources) { return true; }
bool psphotVisualShowPSFStars (pmConfig *config, const pmFPAview *view, pmPSF *psf, psArray *sources) { return true; }
bool psphotVisualShowSatStars (pmConfig *config, const pmFPAview *view, pmPSF *psf, psArray *sources) { return true; }
bool psphotVisualShowPSFModel (pmConfig *config, pmReadout *readout, pmPSF *psf) { return true; }
bool psphotVisualShowFlags (pmConfig *config, const pmFPAview *view, psArray *sources) { return true; }
bool psphotVisualSourceSize (pmConfig *config, const pmFPAview *view, psArray *sources) { return true; }
bool psphotVisualShowResidualImage (pmConfig *config, pmReadout *readout) { return true; }
bool psphotVisualPlotApResid (pmConfig *config, const pmFPAview *view, psArray *sources) { return true; }

# endif
