# include "psphotInternal.h"

// this function displays representative images as the psphot analysis progresses:
// 0 : image, 1 : variance
// 0 : backsub, 1 : variance, 2 : backgnd
// 0 : backsub, 1 : variance, 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->variance, "variance", 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->tmpFlags & PM_SOURCE_TMPF_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
