Index: trunk/psphot/src/Makefile.am =================================================================== --- trunk/psphot/src/Makefile.am (revision 19793) +++ trunk/psphot/src/Makefile.am (revision 19869) @@ -74,4 +74,5 @@ psphotMakeFluxScale.c \ psphotCheckStarDistribution.c \ + psphotVisual.c \ psphotAddNoise.c Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 19793) +++ trunk/psphot/src/psphot.h (revision 19869) @@ -113,5 +113,5 @@ bool psphotRadialPlot (int *kapa, const char *filename, pmSource *source); bool psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe); -bool psphotMosaicSubimage (psImage *outImage, pmSource *source, int Xo, int Yo, int DX, int DY); +bool psphotMosaicSubimage (psImage *outImage, pmSource *source, int Xo, int Yo, int DX, int DY, bool normalize); bool psphotAddWithTest (pmSource *source, bool useState, psMaskType maskVal); @@ -132,4 +132,18 @@ bool psphotAnnuli (pmSource *source, psMetadata *recipe, psMaskType maskVal); bool psphotKron (pmSource *source, psMetadata *recipe, psMaskType maskVal); + +// psphotVisual functions +bool psphotSetVisual (bool mode); +bool psphotVisualShowImage (pmConfig *config, pmReadout *readout); +bool psphotVisualShowBackground (pmConfig *config, const pmFPAview *view, pmReadout *readout); +bool psphotVisualShowSignificance (psImage *image); +bool psphotVisualShowPeaks (pmConfig *config, const pmFPAview *view, pmDetections *detections); +bool psphotVisualShowFootprints (pmConfig *config, const pmFPAview *view, pmDetections *detections); +bool psphotVisualShowMoments (pmConfig *config, const pmFPAview *view, psArray *sources); +bool psphotVisualShowRoughClass (pmConfig *config, const pmFPAview *view, psArray *sources); +bool psphotVisualShowPSFStars (pmConfig *config, const pmFPAview *view, pmPSF *psf, psArray *sources); +bool psphotVisualShowSatStars (pmConfig *config, const pmFPAview *view, pmPSF *psf, psArray *sources); +bool psphotVisualShowLinearFit (pmConfig *config, pmReadout *readout); +bool psphotVisualShowPSFModel (pmConfig *config, pmReadout *readout, pmPSF *psf); // structures & functions to support psf-convolved model fitting Index: trunk/psphot/src/psphotArguments.c =================================================================== --- trunk/psphot/src/psphotArguments.c (revision 19793) +++ trunk/psphot/src/psphotArguments.c (revision 19869) @@ -67,4 +67,10 @@ } + // visual : interactive display mode + if ((N = psArgumentGet (argc, argv, "-visual"))) { + psArgumentRemove (N, &argc, argv); + psphotSetVisual (true); + } + // break : used from recipe throughout psphotReadout if ((N = psArgumentGet (argc, argv, "-break"))) { Index: trunk/psphot/src/psphotBlendFit.c =================================================================== --- trunk/psphot/src/psphotBlendFit.c (revision 19793) +++ trunk/psphot/src/psphotBlendFit.c (revision 19869) @@ -83,4 +83,9 @@ Nfit ++; + // XXX TEST + // if ((fabs(source->peak->x - 1202) < 2) && (fabs(source->peak->y - 1065) < 2)) { + // psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 6); + // } + // try fitting PSFs, then try extended sources // these functions subtract the resulting fitted source (XXX and update the modelFlux?) Index: trunk/psphot/src/psphotFindDetections.c =================================================================== --- trunk/psphot/src/psphotFindDetections.c (revision 19793) +++ trunk/psphot/src/psphotFindDetections.c (revision 19869) @@ -41,4 +41,7 @@ psImage *significance = psphotSignificanceImage (readout, recipe, pass, maskVal); + // display the backsub and backgnd images + psphotVisualShowSignificance (significance); + // detect the peaks in the significance image detections->peaks = psphotFindPeaks (significance, readout, recipe, threshold, NMAX); Index: trunk/psphot/src/psphotFindFootprints.c =================================================================== --- trunk/psphot/src/psphotFindFootprints.c (revision 19793) +++ trunk/psphot/src/psphotFindFootprints.c (revision 19869) @@ -33,4 +33,5 @@ // find the raw footprints & assign the peaks to those footprints psArray *footprints = pmFootprintsFind (significance, threshold, npixMin); + pmFootprintsAssignPeaks(footprints, detections->peaks); // XXX handle the error conditions here Index: trunk/psphot/src/psphotFitSourcesLinear.c =================================================================== --- trunk/psphot/src/psphotFitSourcesLinear.c (revision 19793) +++ trunk/psphot/src/psphotFitSourcesLinear.c (revision 19869) @@ -66,6 +66,8 @@ if (source->type == PM_SOURCE_TYPE_DEFECT) continue; if (source->type == PM_SOURCE_TYPE_SATURATED) continue; - if (source->type == PM_SOURCE_TYPE_STAR && - source->mode & PM_SOURCE_MODE_SATSTAR) continue; + + // if (source->type == PM_SOURCE_TYPE_STAR && + // source->mode & PM_SOURCE_MODE_SATSTAR) continue; + if (final) { if (source->mode & PM_SOURCE_MODE_SUBTRACTED) continue; Index: trunk/psphot/src/psphotMosaicSubimage.c =================================================================== --- trunk/psphot/src/psphotMosaicSubimage.c (revision 19793) +++ trunk/psphot/src/psphotMosaicSubimage.c (revision 19869) @@ -2,5 +2,5 @@ // insert the source image into the outimage at Xo, Yo -bool psphotMosaicSubimage (psImage *outImage, pmSource *source, int Xo, int Yo, int DX, int DY) { +bool psphotMosaicSubimage (psImage *outImage, pmSource *source, int Xo, int Yo, int DX, int DY, bool normalize) { psRegion inRegion, outRegion; @@ -35,5 +35,9 @@ for (int iy = 0; iy < newImage->numRows; iy++) { for (int ix = 0; ix < newImage->numCols; ix++) { - newImage->data.F32[iy][ix] = subImage->data.F32[iy][ix] / peak; + if (normalize) { + newImage->data.F32[iy][ix] = subImage->data.F32[iy][ix] / peak; + } else { + newImage->data.F32[iy][ix] = subImage->data.F32[iy][ix] / peak; + } } } Index: trunk/psphot/src/psphotReadout.c =================================================================== --- trunk/psphot/src/psphotReadout.c (revision 19793) +++ trunk/psphot/src/psphotReadout.c (revision 19869) @@ -36,4 +36,7 @@ } + // set up the Visual display, display the image, weight, mask (ch 1,2,3) + psphotVisualShowImage (config, readout); + // generate a background model (median, smoothed image) if (!psphotModelBackground (config, view, "PSPHOT.INPUT")) { @@ -46,4 +49,7 @@ return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); } + + // display the backsub and backgnd images + psphotVisualShowBackground (config, view, readout); // run a single-model test if desired (exits from here if test is run) @@ -63,4 +69,8 @@ } + // display the backsub and backgnd images + psphotVisualShowPeaks (config, view, detections); + psphotVisualShowFootprints (config, view, detections); + // XXX test write out the footprint image if (0) { @@ -68,11 +78,4 @@ psphotSaveImage (NULL, footprintImage, "footprints.1.fits"); psFree (footprintImage); - } - - // XXX test write out the footprint image - if (0) { - psImage *footprintImage = pmSetFootprintArrayIDs (detections->footprints, true); - psphotSaveImage (NULL, footprintImage, "footprints.1.fits"); - psFree (footprintImage); } @@ -84,4 +87,6 @@ } + psphotVisualShowMoments (config, view, sources); + // find blended neighbors of very saturated stars // XXX merge this with Basic Deblend? @@ -102,4 +107,6 @@ return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); } + + psphotVisualShowRoughClass (config, view, sources); // if we were not supplied a PSF, choose one here @@ -119,4 +126,6 @@ } + psphotVisualShowPSFModel (config, readout, psf); + // Define source fitting parameters for everything that follows PSF fits // This allows different parameters to be used from the PSF fitting @@ -158,4 +167,12 @@ psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); + // XXX we have to place this here because the models are not realized until psphotGuessModels + // or fitted until psphotFitSourcesLinear. There is an example of generating a residual image + // in pmPSFtry.c + + psphotVisualShowPSFStars (config, view, psf, sources); + psphotVisualShowSatStars (config, view, psf, sources); + psphotVisualShowLinearFit (config, readout); + if (0) { FILE *out = fopen ("out.pass1.dat", "w"); @@ -184,6 +201,10 @@ } + // XXX psphotVisualShowSourceSize (config, view, sources); + // non-linear PSF and EXT fit to brighter sources psphotBlendFit (readout, sources, recipe, psf); + + // XXX psphotVisualShowNonLinearFit (config, view, sources); // replace all sources Index: trunk/psphot/src/psphotReplaceUnfit.c =================================================================== --- trunk/psphot/src/psphotReplaceUnfit.c (revision 19793) +++ trunk/psphot/src/psphotReplaceUnfit.c (revision 19869) @@ -71,5 +71,5 @@ } -// add source, if the source has been subtracted (or if we ignore the state) +// add source, if the source has been subtracted; do not modify state bool psphotAddWithTest (pmSource *source, bool useState, psMaskType maskVal) { @@ -78,12 +78,9 @@ if (state && useState) return true; - // replace the model if 1) state says it is missing or 2) useState is false (just do it) - if (!state || !useState) { - pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - } + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); return true; } -// sub source, if the source has been added (or if we ignore the state) +// sub source, if the source has been added; do not modify state bool psphotSubWithTest (pmSource *source, bool useState, psMaskType maskVal) { @@ -92,12 +89,9 @@ if (state && useState) return true; - // replace the model if 1) state says it is missing or 2) useState is false (just do it) - if (!state || !useState) { - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - } + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); return true; } -// add or sub source replace or if the source has +// add or sub source to match recorded state: supply current state as true (add) or false (sub) bool psphotSetState (pmSource *source, bool curState, psMaskType maskVal) { Index: trunk/psphot/src/psphotSavePSFStars.c =================================================================== --- trunk/psphot/src/psphotSavePSFStars.c (revision 19793) +++ trunk/psphot/src/psphotSavePSFStars.c (revision 19869) @@ -100,8 +100,8 @@ psphotAddWithTest (source, true); // replace source if subtracted psphotRadialPlot (&kapa, "radial.plots.ps", source); - psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY); + psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); psphotSubWithTest (source, false); // remove source (force) - psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY); + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); psphotSetState (source, false); // replace source (has been subtracted) @@ -115,8 +115,8 @@ psphotAddWithTest (source, true); // replace source if subtracted psphotRadialPlot (&kapa, "radial.plots.ps", source); - psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY); + psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); psphotSubWithTest (source, false); // remove source (force) - psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY); + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); psphotSetState (source, false); // replace source (has been subtracted) @@ -128,8 +128,8 @@ psphotAddWithTest (source, true); // replace source if subtracted psphotRadialPlot (&kapa, "radial.plots.ps", source); - psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY); + psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); psphotSubWithTest (source, false); // remove source (force) - psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY); + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); psphotSetState (source, false); // replace source (has been subtracted) Index: trunk/psphot/src/psphotSourceFits.c =================================================================== --- trunk/psphot/src/psphotSourceFits.c (revision 19793) +++ trunk/psphot/src/psphotSourceFits.c (revision 19869) @@ -16,5 +16,5 @@ bool psphotFitSummary () { - fprintf (stderr, "fitted %5d psf, %5d blend, %5d ext, %5d dbl : %6.2f sec\n", + psLogMsg ("psphot.pspsf", PS_LOG_INFO, "fitted %5d psf, %5d blend, %5d ext, %5d dbl : %6.2f sec\n", NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits")); return true; @@ -47,4 +47,6 @@ psArray *sourceSet = psArrayAllocEmpty (source->blends->n + 1); psArrayAdd (sourceSet, 16, source); + + psTrace ("psphot", 4, "fitting blended source at %f %f : %f\n", source->peak->xf, source->peak->yf, source->peak->flux); // we need to include all blends in the fit (unless primary is saturated?) @@ -77,4 +79,6 @@ psArrayAdd (modelSet, 16, model); psArrayAdd (sourceSet, 16, blend); + + psTrace ("psphot", 5, "adding source at %f %f : %f\n", blend->peak->xf, blend->peak->yf, blend->peak->flux); // free to avoid double counting model Index: trunk/psphot/src/psphotSourcePlots.c =================================================================== --- trunk/psphot/src/psphotSourcePlots.c (revision 19793) +++ trunk/psphot/src/psphotSourcePlots.c (revision 19869) @@ -113,8 +113,8 @@ psphotAddWithTest (source, true, maskVal); // replace source if subtracted psphotRadialPlot (&kapa, "radial.plots.ps", source); - psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY); + psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); psphotSubWithTest (source, false, maskVal); // remove source (force) - psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY); + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); psphotSetState (source, false, maskVal); // replace source (has been subtracted) @@ -128,8 +128,8 @@ psphotAddWithTest (source, true, maskVal); // replace source if subtracted psphotRadialPlot (&kapa, "radial.plots.ps", source); - psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY); + psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); psphotSubWithTest (source, false, maskVal); // remove source (force) - psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY); + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); psphotSetState (source, false, maskVal); // replace source (has been subtracted) @@ -141,8 +141,8 @@ psphotAddWithTest (source, true, maskVal); // replace source if subtracted psphotRadialPlot (&kapa, "radial.plots.ps", source); - psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY); + psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); psphotSubWithTest (source, false, maskVal); // remove source (force) - psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY); + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); psphotSetState (source, false, maskVal); // replace source (has been subtracted) Index: trunk/psphot/src/psphotVisual.c =================================================================== --- trunk/psphot/src/psphotVisual.c (revision 19869) +++ trunk/psphot/src/psphotVisual.c (revision 19869) @@ -0,0 +1,1092 @@ +# include "psphotInternal.h" +bool psphotVisualPlotMoments (pmConfig *config, const pmFPAview *view, psArray *sources); + +// this function displays representative images as the psphot analysis progresses: +// 0 : image, 1 : weight +// 0 : backsub, 1 : weight, 2 : backgnd +// 0 : backsub, 1 : weight, 2 : signif +// (overlay peaks on images) +// (overlay footprints on images) +// (overlay moments on images) +// (overlay rough class on images) +// 0 : backsub, 1 : psfpos, 2: psfsub +// 0 : backsub, 1 : lin_resid, 2: psfsub + +# if (HAVE_KAPA) +# include + +// 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 psphotVisualScaleImage (int kapaFD, psImage *inImage, const char *name, int channel) { + + KiiImage image; + KapaImageData data; + Coords coords; + + strcpy (coords.ctype, "RA---TAN"); + + psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); + if (!psImageBackground(stats, NULL, inImage, NULL, 0, rng)) { + fprintf (stderr, "failed to get background values\n"); + return false; + } + + image.data2d = inImage->data.F32; + image.Nx = inImage->numCols; + image.Ny = inImage->numRows; + + strcpy (data.name, name); + strcpy (data.file, name); + data.zero = stats->robustMedian - stats->robustStdev; + data.range = 5*stats->robustStdev; + data.logflux = 0; + + KiiSetChannel (kapaFD, channel); + KiiNewPicture2D (kapaFD, &image, &data, &coords); + + psFree (stats); + psFree (rng); + + return true; +} + +bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max) { + + KiiImage image; + KapaImageData data; + Coords coords; + + strcpy (coords.ctype, "RA---TAN"); + + image.data2d = inImage->data.F32; + image.Nx = inImage->numCols; + image.Ny = inImage->numRows; + + strcpy (data.name, name); + strcpy (data.file, name); + data.zero = min; + data.range = max - min; + data.logflux = 0; + + KiiSetChannel (kapaFD, channel); + KiiNewPicture2D (kapaFD, &image, &data, &coords); + + return true; +} + +bool psphotVisualShowImage (pmConfig *config, pmReadout *readout) { + + if (!isVisual) return true; + + if (kapa == -1) { + kapa = KapaOpenNamedSocket ("kapa", "psphot:images"); + if (kapa == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + isVisual = false; + return false; + } + } + + // psphotVisualShowMask (kapa, readout->mask, "mask", 2); + psphotVisualScaleImage (kapa, readout->weight, "weight", 1); + psphotVisualScaleImage (kapa, readout->image, "image", 0); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + return true; +} + +bool psphotVisualShowBackground (pmConfig *config, const pmFPAview *view, pmReadout *readout) { + + pmReadout *backgnd; + + if (!isVisual) return true; + + if (kapa == -1) { + kapa = KapaOpenNamedSocket ("kapa", "psphot:images"); + if (kapa == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + isVisual = false; + return false; + } + } + + bool status = false; + pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKGND"); + + if (file->mode == PM_FPA_MODE_INTERNAL) { + backgnd = file->readout; + } else { + backgnd = pmFPAviewThisReadout (view, file->fpa); + } + + psphotVisualScaleImage (kapa, backgnd->image, "backgnd", 2); + psphotVisualScaleImage (kapa, readout->image, "backsub", 0); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + return true; +} + +bool psphotVisualShowSignificance (psImage *image) { + + if (!isVisual) return true; + + if (kapa == -1) { + kapa = KapaOpenNamedSocket ("kapa", "psphot:images"); + if (kapa == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + isVisual = false; + return false; + } + } + + // XXX test: image->data.F32[10][10] = 10000; + psphotVisualRangeImage (kapa, image, "signif", 2, -1.0, 25.0*25.0); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + return true; +} + +bool psphotVisualShowPeaks (pmConfig *config, const pmFPAview *view, pmDetections *detections) { + + int Noverlay; + KiiOverlay *overlay; + + if (!isVisual) return true; + + if (kapa == -1) { + fprintf (stderr, "kapa not opened, skipping\n"); + return false; + } + + psArray *peaks = detections->peaks; + + // note: this uses the Ohana allocation tools: + // ALLOCATE (overlay, KiiOverlay, 3*peaks->n + 1); + ALLOCATE (overlay, KiiOverlay, peaks->n); + + Noverlay = 0; + for (int i = 0; i < peaks->n; i++) { + + pmPeak *peak = peaks->data[i]; + if (peak == NULL) continue; + + overlay[Noverlay].type = KII_OVERLAY_BOX; + overlay[Noverlay].x = peak->xf; + overlay[Noverlay].y = peak->yf; + overlay[Noverlay].dx = 2.0; + overlay[Noverlay].dy = 2.0; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; + +# if (0) + overlay[Noverlay].type = KII_OVERLAY_BOX; + overlay[Noverlay].x = peak->x; + overlay[Noverlay].y = peak->y; + overlay[Noverlay].dx = 1.0; + overlay[Noverlay].dy = 1.0; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; + + overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + overlay[Noverlay].x = peak->xf; + overlay[Noverlay].y = peak->yf; + overlay[Noverlay].dx = 2.0; + overlay[Noverlay].dy = 2.0; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; +# endif + } + +# if (0) + overlay[Noverlay].type = KII_OVERLAY_BOX; + overlay[Noverlay].x = 10.0; + overlay[Noverlay].y = 10.0; + overlay[Noverlay].dx = 0.5; + overlay[Noverlay].dy = 0.5; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; +# endif + + KiiLoadOverlay (kapa, overlay, Noverlay, "red"); + FREE (overlay); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + return true; +} + +// XXX this is just wrong: outlines entire image +bool psphotVisualShowFootprints (pmConfig *config, const pmFPAview *view, pmDetections *detections) { + + int Noverlay; + KiiOverlay *overlay; + + if (!isVisual) return true; + + if (kapa == -1) { + fprintf (stderr, "kapa not opened, skipping\n"); + return false; + } + + psArray *footprints = detections->footprints; + if (!footprints) return true; + + // note: this uses the Ohana allocation tools: + int NOVERLAY = footprints->n; + ALLOCATE (overlay, KiiOverlay, NOVERLAY); + + Noverlay = 0; + for (int i = 0; i < footprints->n; i++) { + + pmSpan *span = NULL; + + pmFootprint *footprint = footprints->data[i]; + if (footprint == NULL) continue; + if (footprint->spans == NULL) continue; + if (footprint->spans->n < 1) continue; + + // draw the top + span = footprint->spans->data[0]; + overlay[Noverlay].type = KII_OVERLAY_LINE; + overlay[Noverlay].x = span->x0; + overlay[Noverlay].y = span->y; + overlay[Noverlay].dx = span->x1 - span->x0; + overlay[Noverlay].dy = 0; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + + int ys = span->y; + int x0s = span->x0; + int x1s = span->x1; + + // draw the outer span edges + for (int j = 1; j < footprint->spans->n; j++) { + pmSpan *span1 = footprint->spans->data[j]; + + int ye = span1->y; + int x0e = span1->x0; + int x1e = span1->x1; + + // we cannot have two discontinuous spans on the top or bottom, right? (no, probably not right) + // find all of the spans in this row and generate x0e, x01: + for (int k = j + 1; k < footprint->spans->n; k++) { + pmSpan *span2 = footprint->spans->data[k]; + if (span2->y > span1->y) break; + x0e = PS_MIN (x0e, span2->x0); + x1e = PS_MAX (x1e, span2->x1); + j++; + } + + overlay[Noverlay].type = KII_OVERLAY_LINE; + overlay[Noverlay].x = x0s; + overlay[Noverlay].y = ys; + overlay[Noverlay].dx = x0e - x0s; + overlay[Noverlay].dy = ye - ys; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + + overlay[Noverlay].type = KII_OVERLAY_LINE; + overlay[Noverlay].x = x1s; + overlay[Noverlay].y = ys; + overlay[Noverlay].dx = x1e - x1s; + overlay[Noverlay].dy = ye - ys; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + + ys = ye; + x0s = x0e; + x1s = x1e; + } + + // draw the bottom + span = footprint->spans->data[footprint->spans->n - 1]; + overlay[Noverlay].type = KII_OVERLAY_LINE; + overlay[Noverlay].x = span->x0; + overlay[Noverlay].y = span->y; + overlay[Noverlay].dx = span->x1 - span->x0; + overlay[Noverlay].dy = 0; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + } + + KiiLoadOverlay (kapa, overlay, Noverlay, "blue"); + FREE (overlay); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + return true; +} + +bool psphotVisualShowMoments (pmConfig *config, const pmFPAview *view, psArray *sources) { + + int Noverlay; + KiiOverlay *overlay; + + psEllipseMoments emoments; + psEllipseAxes axes; + + if (!isVisual) return true; + + if (kapa == -1) { + fprintf (stderr, "kapa not opened, skipping\n"); + return false; + } + + // note: this uses the Ohana allocation tools: + ALLOCATE (overlay, KiiOverlay, sources->n); + + Noverlay = 0; + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + if (source == NULL) continue; + + pmMoments *moments = source->moments; + if (moments == NULL) continue; + + overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + overlay[Noverlay].x = moments->x; + overlay[Noverlay].y = moments->y; + + emoments.x2 = moments->Sx; + emoments.y2 = moments->Sy; + emoments.xy = moments->Sxy; + + axes = psEllipseMomentsToAxes (emoments, 20.0); + + overlay[Noverlay].dx = 2.0*axes.major; + overlay[Noverlay].dy = 2.0*axes.minor; + overlay[Noverlay].angle = axes.theta * PS_DEG_RAD; + overlay[Noverlay].text = NULL; + Noverlay ++; + } + + KiiLoadOverlay (kapa, overlay, Noverlay, "yellow"); + FREE (overlay); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + + return true; +} + +bool psphotVisualPlotMoments (pmConfig *config, const pmFPAview *view, psArray *sources) { + + bool status; + Graphdata graphdata; + + if (!isVisual) return true; + + if (kapa3 == -1) { + kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots"); + if (kapa3 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + isVisual = false; + return false; + } + } + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); + if (!recipe) { + fprintf (stderr, "missing recipe, skipping moments plot\n"); + return false; + } + + KapaClearPlots (kapa3); + KapaInitGraph (&graphdata); + + float psfX = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.X"); + float psfY = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.Y"); + float psfdX = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DX"); + float psfdY = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DY"); + float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA"); + float Rx = psfdX * PSF_CLUMP_NSIGMA; + float Ry = psfdY * PSF_CLUMP_NSIGMA; + + // examine sources to set data range + graphdata.xmin = -0.05; + graphdata.ymin = -0.05; + graphdata.xmax = 2.0*psfX; + graphdata.ymax = 2.0*psfY; + KapaSetLimits (kapa3, &graphdata); + + KapaSetFont (kapa3, "helvetica", 14); + KapaBox (kapa3, &graphdata); + KapaSendLabel (kapa3, "&ss&h_x| (pixels)", KAPA_LABEL_XM); + KapaSendLabel (kapa3, "&ss&h_y| (pixels)", KAPA_LABEL_YM); + + psVector *xBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *yBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *xFaint = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *yFaint = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + + // construct the vectors + int nB = 0; + int nF = 0; + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (source->moments == NULL) + continue; + + xFaint->data.F32[nF] = source->moments->Sx; + yFaint->data.F32[nF] = source->moments->Sy; + nF++; + + // XXX make this a user-defined cutoff + if (source->moments->SN < 50) + continue; + + xBright->data.F32[nB] = source->moments->Sx; + yBright->data.F32[nB] = source->moments->Sy; + nB++; + } + xFaint->n = nF; + yFaint->n = nF; + + xBright->n = nB; + yBright->n = nB; + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 0; + graphdata.size = 0.3; + graphdata.style = 2; + KapaPrepPlot (kapa3, nF, &graphdata); + KapaPlotVector (kapa3, nF, xFaint->data.F32, "x"); + KapaPlotVector (kapa3, nF, yFaint->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 2; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (kapa3, nB, &graphdata); + KapaPlotVector (kapa3, nB, xBright->data.F32, "x"); + KapaPlotVector (kapa3, nB, yBright->data.F32, "y"); + + // draw a circle centered on psfX,Y with size of the psf limit + psVector *xLimit = psVectorAlloc (120, PS_TYPE_F32); + psVector *yLimit = psVectorAlloc (120, PS_TYPE_F32); + for (int i = 0; i < xLimit->n; i++) { + xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX; + yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY; + } + graphdata.color = KapaColorByName ("blue"); + graphdata.style = 0; + KapaPrepPlot (kapa3, xLimit->n, &graphdata); + KapaPlotVector (kapa3, xLimit->n, xLimit->data.F32, "x"); + KapaPlotVector (kapa3, yLimit->n, yLimit->data.F32, "y"); + psFree (xLimit); + psFree (yLimit); + +# if (0) + // *** make a histogram of the source counts in the x and y directions + psHistogram *nX = psHistogramAlloc (graphdata.xmin, graphdata.xmax, 50.0); + psHistogram *nY = psHistogramAlloc (graphdata.ymin, graphdata.ymax, 50.0); + psVectorHistogram (nX, xFaint, NULL, NULL, 0); + psVectorHistogram (nY, yFaint, NULL, NULL, 0); + psVector *dX = psVectorAlloc (nX->nums->n, PS_TYPE_F32); + psVector *vX = psVectorAlloc (nX->nums->n, PS_TYPE_F32); + psVector *dY = psVectorAlloc (nY->nums->n, PS_TYPE_F32); + psVector *vY = psVectorAlloc (nY->nums->n, PS_TYPE_F32); + for (int i = 0; i < nX->nums->n; i++) { + dX->data.F32[i] = nX->nums->data.S32[i]; + vX->data.F32[i] = 0.5*(nX->bounds->data.F32[i] + nX->bounds->data.F32[i+1]); + } + for (int i = 0; i < nY->nums->n; i++) { + dY->data.F32[i] = nY->nums->data.S32[i]; + vY->data.F32[i] = 0.5*(nY->bounds->data.F32[i] + nY->bounds->data.F32[i+1]); + } + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 0; + graphdata.size = 0.0; + graphdata.style = 0; + KapaPrepPlot (kapa3, dX->n, &graphdata); + KapaPlotVector (kapa3, dX->n, dX->data.F32, "x"); + KapaPlotVector (kapa3, vX->n, vX->data.F32, "y"); + + psFree (nX); + psFree (dX); + psFree (vX); + + psFree (nY); + psFree (dY); + psFree (vY); +# endif + + psFree (xBright); + psFree (yBright); + psFree (xFaint); + psFree (yFaint); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + return true; +} + +// assumes 'kapa' value is checked and set +bool psphotVisualShowRoughClass_Single (pmConfig *config, const pmFPAview *view, psArray *sources, pmSourceType type, pmSourceMode mode, char *color) { + + int Noverlay; + KiiOverlay *overlay; + + psEllipseMoments emoments; + psEllipseAxes axes; + + // note: this uses the Ohana allocation tools: + ALLOCATE (overlay, KiiOverlay, sources->n); + + Noverlay = 0; + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + if (source == NULL) continue; + + if (source->type != type) continue; + if (mode && !(source->mode & mode)) continue; + + pmMoments *moments = source->moments; + if (moments == NULL) continue; + + overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + overlay[Noverlay].x = moments->x; + overlay[Noverlay].y = moments->y; + + emoments.x2 = moments->Sx; + emoments.y2 = moments->Sy; + emoments.xy = moments->Sxy; + + axes = psEllipseMomentsToAxes (emoments, 20.0); + + overlay[Noverlay].dx = 2.0*axes.major; + overlay[Noverlay].dy = 2.0*axes.minor; + overlay[Noverlay].angle = axes.theta * PS_DEG_RAD; + overlay[Noverlay].text = NULL; + Noverlay ++; + } + + KiiLoadOverlay (kapa, overlay, Noverlay, color); + FREE (overlay); + + return true; +} + +bool psphotVisualShowRoughClass (pmConfig *config, const pmFPAview *view, psArray *sources) { + + psphotVisualPlotMoments (config, view, sources); + + if (!isVisual) return true; + + if (kapa == -1) { + fprintf (stderr, "kapa not opened, skipping\n"); + return false; + } + + KiiEraseOverlay (kapa, "yellow"); // moments + + psphotVisualShowRoughClass_Single (config, view, sources, PM_SOURCE_TYPE_STAR, 0, "red"); + psphotVisualShowRoughClass_Single (config, view, sources, PM_SOURCE_TYPE_EXTENDED, 0, "blue"); + psphotVisualShowRoughClass_Single (config, view, sources, PM_SOURCE_TYPE_DEFECT, 0, "blue"); + psphotVisualShowRoughClass_Single (config, view, sources, PM_SOURCE_TYPE_SATURATED, 0, "red"); + psphotVisualShowRoughClass_Single (config, view, sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_PSFSTAR, "yellow"); + psphotVisualShowRoughClass_Single (config, view, sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_SATSTAR, "green"); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "red: STAR or SAT AREA; blue: EXTENDED or DEFECT; green: SATSTAR; yellow: PSFSTAR\n"); + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + return true; +} + +bool psphotVisualShowPSFStars (pmConfig *config, const pmFPAview *view, pmPSF *psf, psArray *sources) { + + bool status; + + if (!isVisual) return true; + + if (kapa2 == -1) { + kapa2 = KapaOpenNamedSocket ("kapa", "psphot:psfstars"); + if (kapa2 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + isVisual = false; + return false; + } + } + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); + if (!recipe) { + psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE); + return false; + } + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // the source images are written to an image 10x the size of a PSF object + // float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS"); + // PS_ASSERT (status, false); + + int DX = 21; + int DY = 21; + + // examine PSF sources in S/N order (brightest first) + sources = psArraySort (sources, pmSourceSortBySN); + + // counters to track the size of the image and area used in a row + int dX = 0; // starting corner of next box + int dY = 0; // height of row so far + int NX = 20*DX; // full width of output image + int NY = 0; // total height of output image + + // first, examine the PSF stars: + // - determine bounding boxes for summary image + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + bool keep = false; + keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR); + if (!keep) continue; + + // how does this subimage get placed into the output image? + // DX = source->pixels->numCols + // DY = source->pixels->numRows + + if (dX + DX > NX) { + // too wide for the rest of this row + if (dX == 0) { + // alone on this row + NY += DY; + dX = 0; + dY = 0; + } else { + // start the next row + NY += dY; + dX = DX; + dY = DY; + } + } else { + // extend this row + dX += DX; + dY = PS_MAX (dY, DY); + } + } + NY += DY; + + // allocate output image + psImage *outpos = psImageAlloc (NX, NY, PS_TYPE_F32); + psImage *outsub = psImageAlloc (NX, NY, PS_TYPE_F32); + + int Xo = 0; // starting corner of next box + int Yo = 0; // starting corner of next box + dY = 0; // height of row so far + + int nPSF = 0; + + // next, examine the PSF stars: + // - create output image array + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + bool keep = false; + if (source->mode & PM_SOURCE_MODE_PSFSTAR) { + nPSF ++; + keep = true; + } + if (!keep) continue; + + if (Xo + DX > NX) { + // too wide for the rest of this row + if (Xo == 0) { + // place source alone on this row + psphotAddWithTest (source, true, maskVal); // replace source if subtracted + psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); + + psphotSubWithTest (source, false, maskVal); // remove source (force) + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); + psphotSetState (source, false, maskVal); // reset source Add/Sub state to recorded + + Yo += DY; + Xo = 0; + dY = 0; + } else { + // start the next row + Yo += dY; + Xo = 0; + + psphotAddWithTest (source, true, maskVal); // replace source if subtracted + psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); + + psphotSubWithTest (source, false, maskVal); // remove source (force) + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); + psphotSetState (source, false, maskVal); // replace source (has been subtracted) + + Xo = DX; + dY = DY; + } + } else { + // extend this row + psphotAddWithTest (source, true, maskVal); // replace source if subtracted + psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); + + psphotSubWithTest (source, false, maskVal); // remove source (force) + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); + psphotSetState (source, false, maskVal); // replace source (has been subtracted) + + Xo += DX; + dY = PS_MAX (dY, DY); + } + } + + psphotVisualScaleImage (kapa2, outpos, "psfpos", 0); + psphotVisualScaleImage (kapa2, outsub, "psfsub", 1); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + + psFree (outpos); + psFree (outsub); + return true; +} + +bool psphotVisualShowSatStars (pmConfig *config, const pmFPAview *view, pmPSF *psf, psArray *sources) { + + bool status; + + if (!isVisual) return true; + + if (kapa2 == -1) { + kapa2 = KapaOpenNamedSocket ("kapa", "psphot:images"); + if (kapa2 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + isVisual = false; + return false; + } + } + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); + if (!recipe) { + psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE); + return false; + } + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // the source images are written to an image 10x the size of a PSF object + // float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS"); + // PS_ASSERT (status, false); + + int DX = 41; + int DY = 41; + + // examine PSF sources in S/N order (brightest first) + sources = psArraySort (sources, pmSourceSortBySN); + + // counters to track the size of the image and area used in a row + int dX = 0; // starting corner of next box + int dY = 0; // height of row so far + int NX = 10*DX; // full width of output image + int NY = 0; // total height of output image + + // first, examine the PSF and SAT stars: + // - determine bounding boxes for summary image + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + bool keep = false; + keep |= (source->mode & PM_SOURCE_MODE_SATSTAR); + if (!keep) continue; + + // how does this subimage get placed into the output image? + // DX = source->pixels->numCols + // DY = source->pixels->numRows + + if (dX + DX > NX) { + // too wide for the rest of this row + if (dX == 0) { + // alone on this row + NY += DY; + dX = 0; + dY = 0; + } else { + // start the next row + NY += dY; + dX = DX; + dY = DY; + } + } else { + // extend this row + dX += DX; + dY = PS_MAX (dY, DY); + } + } + NY += DY; + + // allocate output image + psImage *outsat = psImageAlloc (NX, NY, PS_TYPE_F32); + + int Xo = 0; // starting corner of next box + int Yo = 0; // starting corner of next box + dY = 0; // height of row so far + + int nSAT = 0; + + // next, examine the SAT stars: + // - create output image array + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + bool keep = false; + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + nSAT ++; + keep = true; + } + if (!keep) continue; + + if (Xo + DX > NX) { + // too wide for the rest of this row + if (Xo == 0) { + // place source alone on this row + psphotAddWithTest (source, true, maskVal); // replace source if subtracted + psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); + psphotSetState (source, true, maskVal); // reset source Add/Sub state to recorded + + Yo += DY; + Xo = 0; + dY = 0; + } else { + // start the next row + Yo += dY; + Xo = 0; + psphotAddWithTest (source, true, maskVal); // replace source if subtracted + psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); + psphotSetState (source, true, maskVal); // replace source (has been subtracted) + + Xo = DX; + dY = DY; + } + } else { + // extend this row + psphotAddWithTest (source, true, maskVal); // replace source if subtracted + psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); + psphotSetState (source, true, maskVal); // replace source (has been subtracted) + + Xo += DX; + dY = PS_MAX (dY, DY); + } + } + + psphotVisualScaleImage (kapa2, outsat, "satstar", 2); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + + psFree (outsat); + return true; +} + +bool psphotVisualShowPSFModel (pmConfig *config, pmReadout *readout, pmPSF *psf) { + + if (!isVisual) return true; + + if (kapa2 == -1) { + kapa2 = KapaOpenNamedSocket ("kapa", "psphot:psfstars"); + if (kapa2 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + isVisual = false; + return false; + } + } + + int DX = 64; + int DY = 64; + + psImage *psfMosaic = psImageAlloc (3*DX, 3*DY, PS_TYPE_F32); + psImageInit (psfMosaic, 0.0); + + pmModel *modelRef = pmModelAlloc(psf->type); + +# if (1) + // generate a fake model at each of the 3x3 image grid positions + for (int x = -1; x <= +1; x ++) { + for (int y = -1; y <= +1; y ++) { + // use the center of the center pixel of the image + float xc = (0.5 + 0.45*x)*readout->image->numCols + readout->image->col0; + float yc = (0.5 + 0.45*y)*readout->image->numRows + readout->image->row0; + + // assign the x and y coords to the image center + // create an object with center intensity of 1000 + modelRef->params->data.F32[PM_PAR_SKY] = 0; + modelRef->params->data.F32[PM_PAR_I0] = 1000; + modelRef->params->data.F32[PM_PAR_XPOS] = xc; + modelRef->params->data.F32[PM_PAR_YPOS] = yc; + + // create modelPSF from this model + pmModel *model = pmModelFromPSF (modelRef, psf); + + // place the reference object in the image center + // no need to mask the source here + // XXX should we measure this for the analytical model only or the full model? + pmModelAddWithOffset (psfMosaic, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_CENTER, 0, -x*DX, -y*DY); + } + } +# else + { + // use the center of the center pixel of the image + float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5; + float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5; + + // assign the x and y coords to the image center + // create an object with center intensity of 1000 + modelRef->params->data.F32[PM_PAR_SKY] = 0; + modelRef->params->data.F32[PM_PAR_I0] = 1000; + modelRef->params->data.F32[PM_PAR_XPOS] = xc; + modelRef->params->data.F32[PM_PAR_YPOS] = yc; + + // create modelPSF from this model + pmModel *model = pmModelFromPSF (modelRef, psf); + + // place the reference object in the image center + // no need to mask the source here + // XXX should we measure this for the analytical model only or the full model? + pmModelAddWithOffset (psfMosaic, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_CENTER, 0, 32.0, 16.0); + } +# endif + + psImage *psfLogFlux = (psImage *) psUnaryOp (NULL, psfMosaic, "log"); + psphotVisualRangeImage (kapa2, psfLogFlux, "psf_mosaic", 1, -2.0, 3.0); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + return true; +} + +bool psphotVisualShowLinearFit (pmConfig *config, pmReadout *readout) { + + if (!isVisual) return true; + + if (kapa == -1) { + kapa = KapaOpenNamedSocket ("kapa", "psphot:images"); + if (kapa == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + isVisual = false; + return false; + } + } + + psphotVisualScaleImage (kapa, readout->image, "lin_resid", 1); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + return true; +} + +# else + +bool psphotSetVisual (bool mode){} +bool psphotVisualShowImage (pmConfig *config, pmReadout *readout){} +bool psphotVisualShowBackground (pmConfig *config, const pmFPAview *view, pmReadout *readout){} +bool psphotVisualShowSignificance (psImage *image){} +bool psphotVisualShowPeaks (pmConfig *config, const pmFPAview *view, pmDetections *detections){} +bool psphotVisualShowFootprints (pmConfig *config, const pmFPAview *view, pmDetections *detections){} +bool psphotVisualShowMoments (pmConfig *config, const pmFPAview *view, psArray *sources){} +bool psphotVisualShowRoughClass (pmConfig *config, const pmFPAview *view, psArray *sources){} +bool psphotVisualShowPSF (pmConfig *config, const pmFPAview *view, pmPSF *psf, psArray *sources){} +bool psphotVisualShowLinearFit (pmConfig *config, pmReadout *readout){} + +# endif