Index: branches/simmosaic_branches/psphot/src/psphotVisual.c =================================================================== --- branches/simmosaic_branches/psphot/src/psphotVisual.c (revision 24860) +++ branches/simmosaic_branches/psphot/src/psphotVisual.c (revision 27839) @@ -15,11 +15,62 @@ # include +bool pmVisualLimitsFromVectors (Graphdata *graphdata, psVector *xVec, psVector *yVec); + // functions used to visualize the analysis as it goes // these are invoked by the -visual options -static int kapa = -1; +static int kapa1 = -1; static int kapa2 = -1; static int kapa3 = -1; +int psphotKapaChannel (int channel) { + + switch (channel) { + case 1: + if (kapa1 == -1) { + kapa1 = KapaOpenNamedSocket ("kapa", "psphot:images"); + if (kapa1 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + } + } + return kapa1; + case 2: + if (kapa2 == -1) { + kapa2 = KapaOpenNamedSocket ("kapa", "psphot:plots"); + if (kapa2 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + } + } + return kapa2; + case 3: + if (kapa3 == -1) { + kapa3 = KapaOpenNamedSocket ("kapa", "psphot:stamps"); + if (kapa3 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + } + } + return kapa3; + default: + psAbort ("unknown kapa channel"); + } + psAbort ("unknown kapa channel"); +} + +bool psphotVisualEraseOverlays (int channel, char *overlay) { + + int myKapa = psphotKapaChannel (channel); + if (!(strcasecmp (overlay, "all"))) { + KiiEraseOverlay (myKapa, "red"); + KiiEraseOverlay (myKapa, "green"); + KiiEraseOverlay (myKapa, "blue"); + KiiEraseOverlay (myKapa, "yellow"); + return true; + } + KiiEraseOverlay (myKapa, overlay); + return true; +} bool psphotVisualShowMask (int kapaFD, psImage *inImage, const char *name, int channel) { @@ -69,5 +120,5 @@ } -bool psphotVisualScaleImage (int kapaFD, psImage *inImage, const char *name, int channel) { +bool psphotVisualScaleImage (int kapaFD, psImage *inImage, psImage *inMask, const char *name, int channel) { KiiImage image; @@ -79,5 +130,5 @@ psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); - if (!psImageBackground(stats, NULL, inImage, NULL, 0, rng)) { + if (!psImageBackground(stats, NULL, inImage, inMask, 0xffff, rng)) { fprintf (stderr, "failed to get background values\n"); return false; @@ -131,24 +182,12 @@ if (!pmVisualIsVisual()) return true; - if (kapa == -1) { - kapa = KapaOpenNamedSocket ("kapa", "psphot:images"); - if (kapa == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(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"); - } + int kapa = psphotKapaChannel (1); + if (kapa == -1) return false; + + psphotVisualShowMask (kapa, readout->mask, "mask", 2); + psphotVisualScaleImage (kapa, readout->variance, readout->mask, "variance", 1); + psphotVisualScaleImage (kapa, readout->image, readout->mask, "image", 0); + + pmVisualAskUser(NULL); return true; } @@ -160,12 +199,6 @@ if (!pmVisualIsVisual()) return true; - if (kapa == -1) { - kapa = KapaOpenNamedSocket ("kapa", "psphot:images"); - if (kapa == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } + int kapa = psphotKapaChannel (1); + if (kapa == -1) return false; bool status = false; @@ -178,40 +211,21 @@ } - 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) { + psphotVisualScaleImage (kapa, backgnd->image, readout->mask, "backgnd", 2); + psphotVisualScaleImage (kapa, readout->image, readout->mask, "backsub", 0); + + pmVisualAskUser(NULL); + return true; +} + +bool psphotVisualShowSignificance (psImage *image, float min, float max) { if (!pmVisualIsVisual()) return true; - if (kapa == -1) { - kapa = KapaOpenNamedSocket ("kapa", "psphot:images"); - if (kapa == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(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"); - } + int kapa = psphotKapaChannel (1); + if (kapa == -1) return false; + + psphotVisualRangeImage (kapa, image, "signif", 2, min, max); + + pmVisualAskUser(NULL); return true; } @@ -224,8 +238,6 @@ if (!pmVisualIsVisual()) return true; - if (kapa == -1) { - fprintf (stderr, "kapa not opened, skipping\n"); - return false; - } + int kapa = psphotKapaChannel (1); + if (kapa == -1) return false; psArray *peaks = detections->peaks; @@ -233,5 +245,5 @@ // note: this uses the Ohana allocation tools: // ALLOCATE (overlay, KiiOverlay, 3*peaks->n + 1); - ALLOCATE (overlay, KiiOverlay, peaks->n); + ALLOCATE (overlay, KiiOverlay, peaks->n + 2); Noverlay = 0; @@ -249,47 +261,10 @@ 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"); - } + pmVisualAskUser(NULL); return true; } @@ -302,8 +277,6 @@ if (!pmVisualIsVisual()) return true; - if (kapa == -1) { - fprintf (stderr, "kapa not opened, skipping\n"); - return false; - } + int kapa = psphotKapaChannel (1); + if (kapa == -1) return false; psArray *footprints = detections->footprints; @@ -325,4 +298,5 @@ // draw the top + // XXX need to allow top (and bottom) to have more than one span span = footprint->spans->data[0]; overlay[Noverlay].type = KII_OVERLAY_LINE; @@ -399,11 +373,5 @@ 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"); - } + pmVisualAskUser(NULL); return true; } @@ -419,8 +387,9 @@ if (!pmVisualIsVisual()) return true; - if (kapa == -1) { - fprintf (stderr, "kapa not opened, skipping\n"); - return false; - } + int kapa = psphotKapaChannel (1); + if (kapa == -1) return false; + + // XXX mark the different source classes with different color/shape dots + // XXX are moments S/N and peak S/N consistent? // note: this uses the Ohana allocation tools: @@ -448,5 +417,7 @@ 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].angle = axes.theta * PS_DEG_RAD; + overlay[Noverlay].text = NULL; Noverlay ++; @@ -456,16 +427,9 @@ 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) { + pmVisualAskUser(NULL); + return true; +} + +bool psphotVisualPlotMoments (psMetadata *recipe, psMetadata *analysis, psArray *sources) { bool status; @@ -475,46 +439,48 @@ if (!pmVisualIsVisual()) return true; - if (kapa3 == -1) { - kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots"); - if (kapa3 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } - - KapaClearPlots (kapa3); + int myKapa = psphotKapaChannel (2); + if (myKapa == -1) return false; + + KapaClearPlots (myKapa); KapaInitGraph (&graphdata); - KapaSetFont (kapa3, "courier", 14); + KapaSetFont (myKapa, "courier", 14); + + section.bg = KapaColorByName ("none"); // XXX probably should be 'none' float SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); // select the max psfX,Y values for the plot limits - float Xmin = 0.0, Xmax = 0.0; - float Ymin = 0.0, Ymax = 0.0; + float Xmin = 1000.0, Xmax = 0.0; + float Ymin = 1000.0, Ymax = 0.0; { - int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS"); + int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS"); 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"); + psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, 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 X0 = psfX - 4.0*psfdX; - float X1 = psfX + 4.0*psfdX; - float Y0 = psfY - 4.0*psfdY; - float Y1 = psfY + 4.0*psfdY; - - if (isfinite(X0)) { Xmin = PS_MAX(Xmin, X0); } - if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); } - if (isfinite(Y0)) { Ymin = PS_MAX(Ymin, Y0); } - if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); } - } - } + float X0 = psfX - 4.0*psfdX; + float X1 = psfX + 4.0*psfdX; + float Y0 = psfY - 4.0*psfdY; + float Y1 = psfY + 4.0*psfdY; + + if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); } + if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); } + if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); } + if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); } + } + } + Xmin = PS_MAX(Xmin, -0.1); + Ymin = PS_MAX(Ymin, -0.1); + + // XXX test: hardwire plot limits + // Xmin = -0.1; Ymin = -0.1; + // Xmax = 20.1; Ymax = 20.1; // storage vectors for data to be plotted @@ -564,5 +530,5 @@ section.y = 0.00; section.name = psStringCopy ("MxxMyy"); - KapaSetSection (kapa3, §ion); + KapaSetSection (myKapa, §ion); psFree (section.name); @@ -572,9 +538,14 @@ graphdata.xmax = Xmax; graphdata.ymax = Ymax; - KapaSetLimits (kapa3, &graphdata); - - KapaBox (kapa3, &graphdata); - KapaSendLabel (kapa3, "M_xx| (pixels)", KAPA_LABEL_XM); - KapaSendLabel (kapa3, "M_yy| (pixels)", KAPA_LABEL_YM); + KapaSetLimits (myKapa, &graphdata); + + graphdata.padXm = NAN; + graphdata.padYm = NAN; + graphdata.padXp = 0.5; + graphdata.padYp = 0.5; + KapaBox (myKapa, &graphdata); + + KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM); + KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM); graphdata.color = KapaColorByName ("black"); @@ -582,7 +553,8 @@ 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"); + KapaPrepPlot (myKapa, nF, &graphdata); + + KapaPlotVector (myKapa, nF, xFaint->data.F32, "x"); + KapaPlotVector (myKapa, nF, yFaint->data.F32, "y"); graphdata.color = KapaColorByName ("red"); @@ -590,7 +562,7 @@ 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"); + KapaPrepPlot (myKapa, nB, &graphdata); + KapaPlotVector (myKapa, nB, xBright->data.F32, "x"); + KapaPlotVector (myKapa, nB, yBright->data.F32, "y"); // second section: MagMyy @@ -598,7 +570,7 @@ section.dy = 0.25; section.x = 0.00; - section.y = 0.80; + section.y = 0.75; section.name = psStringCopy ("MagMyy"); - KapaSetSection (kapa3, §ion); + KapaSetSection (myKapa, §ion); psFree (section.name); @@ -608,10 +580,14 @@ graphdata.ymin = Ymin; graphdata.ymax = Ymax; - KapaSetLimits (kapa3, &graphdata); - + KapaSetLimits (myKapa, &graphdata); + + graphdata.padXm = 0.5; + graphdata.padYm = NAN; + graphdata.padXp = NAN; + graphdata.padYp = 0.5; strcpy (graphdata.labels, "0210"); - KapaBox (kapa3, &graphdata); - KapaSendLabel (kapa3, "inst mag", KAPA_LABEL_XP); - KapaSendLabel (kapa3, "M_yy| (pixels)", KAPA_LABEL_YM); + KapaBox (myKapa, &graphdata); + KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_XP); + KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM); graphdata.color = KapaColorByName ("black"); @@ -619,7 +595,7 @@ graphdata.size = 0.3; graphdata.style = 2; - KapaPrepPlot (kapa3, nF, &graphdata); - KapaPlotVector (kapa3, nF, mFaint->data.F32, "x"); - KapaPlotVector (kapa3, nF, yFaint->data.F32, "y"); + KapaPrepPlot (myKapa, nF, &graphdata); + KapaPlotVector (myKapa, nF, mFaint->data.F32, "x"); + KapaPlotVector (myKapa, nF, yFaint->data.F32, "y"); graphdata.color = KapaColorByName ("red"); @@ -627,15 +603,15 @@ graphdata.size = 0.5; graphdata.style = 2; - KapaPrepPlot (kapa3, nB, &graphdata); - KapaPlotVector (kapa3, nB, mBright->data.F32, "x"); - KapaPlotVector (kapa3, nB, yBright->data.F32, "y"); + KapaPrepPlot (myKapa, nB, &graphdata); + KapaPlotVector (myKapa, nB, mBright->data.F32, "x"); + KapaPlotVector (myKapa, nB, yBright->data.F32, "y"); // third section: MagMxx section.dx = 0.25; section.dy = 0.75; - section.x = 0.80; + section.x = 0.75; section.y = 0.00; section.name = psStringCopy ("MagMxx"); - KapaSetSection (kapa3, §ion); + KapaSetSection (myKapa, §ion); psFree (section.name); @@ -645,10 +621,14 @@ graphdata.ymin = -7.9; graphdata.ymax = -17.1; - KapaSetLimits (kapa3, &graphdata); - + KapaSetLimits (myKapa, &graphdata); + + graphdata.padXm = NAN; + graphdata.padYm = 0.5; + graphdata.padXp = 0.5; + graphdata.padYp = NAN; strcpy (graphdata.labels, "2001"); - KapaBox (kapa3, &graphdata); - KapaSendLabel (kapa3, "M_xx| (pixels)", KAPA_LABEL_XM); - KapaSendLabel (kapa3, "inst mag", KAPA_LABEL_YP); + KapaBox (myKapa, &graphdata); + KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM); + KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_YP); graphdata.color = KapaColorByName ("black"); @@ -656,7 +636,7 @@ graphdata.size = 0.3; graphdata.style = 2; - KapaPrepPlot (kapa3, nF, &graphdata); - KapaPlotVector (kapa3, nF, xFaint->data.F32, "x"); - KapaPlotVector (kapa3, nF, mFaint->data.F32, "y"); + KapaPrepPlot (myKapa, nF, &graphdata); + KapaPlotVector (myKapa, nF, xFaint->data.F32, "x"); + KapaPlotVector (myKapa, nF, mFaint->data.F32, "y"); graphdata.color = KapaColorByName ("red"); @@ -664,11 +644,11 @@ graphdata.size = 0.5; graphdata.style = 2; - KapaPrepPlot (kapa3, nB, &graphdata); - KapaPlotVector (kapa3, nB, xBright->data.F32, "x"); - KapaPlotVector (kapa3, nB, mBright->data.F32, "y"); + KapaPrepPlot (myKapa, nB, &graphdata); + KapaPlotVector (myKapa, nB, xBright->data.F32, "x"); + KapaPlotVector (myKapa, nB, mBright->data.F32, "y"); // draw N circles to outline the clumps { - KapaSelectSection (kapa3, "MxxMyy"); + KapaSelectSection (myKapa, "MxxMyy"); // draw a circle centered on psfX,Y with size of the psf limit @@ -676,5 +656,5 @@ psVector *yLimit = psVectorAlloc (120, PS_TYPE_F32); - int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS"); + int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS"); float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA"); @@ -682,9 +662,9 @@ graphdata.style = 0; - graphdata.xmin = Xmin; - graphdata.ymin = Ymin; - graphdata.xmax = Xmax; - graphdata.ymax = Ymax; - KapaSetLimits (kapa3, &graphdata); + graphdata.xmin = Xmin; + graphdata.ymin = Ymin; + graphdata.xmax = Xmax; + graphdata.ymax = Ymax; + KapaSetLimits (myKapa, &graphdata); for (int n = 0; n < nRegions; n++) { @@ -692,5 +672,5 @@ char regionName[64]; snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n); - psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName); + psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); @@ -705,47 +685,11 @@ 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"); + KapaPrepPlot (myKapa, xLimit->n, &graphdata); + KapaPlotVector (myKapa, xLimit->n, xLimit->data.F32, "x"); + KapaPlotVector (myKapa, 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); @@ -756,16 +700,10 @@ psFree (mFaint); - // 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"); - } + pmVisualAskUser(NULL); return true; } // assumes 'kapa' value is checked and set -bool psphotVisualShowRoughClass_Single (psArray *sources, pmSourceType type, pmSourceMode mode, char *color) { +bool psphotVisualShowRoughClass_Single (int myKapa, psArray *sources, pmSourceType type, pmSourceMode mode, char *color) { int Noverlay; @@ -802,10 +740,10 @@ overlay[Noverlay].dx = 2.0*axes.major; overlay[Noverlay].dy = 2.0*axes.minor; - overlay[Noverlay].angle = -axes.theta * PS_DEG_RAD; + overlay[Noverlay].angle = axes.theta * PS_DEG_RAD; overlay[Noverlay].text = NULL; Noverlay ++; } - KiiLoadOverlay (kapa, overlay, Noverlay, color); + KiiLoadOverlay (myKapa, overlay, Noverlay, color); FREE (overlay); @@ -817,26 +755,18 @@ if (!pmVisualIsVisual()) 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]; + int myKapa = psphotKapaChannel (1); + if (myKapa == -1) return false; + + KiiEraseOverlay (myKapa, "yellow"); // moments + + psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_STAR, 0, "red"); + psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_EXTENDED, 0, "blue"); + psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_DEFECT, 0, "blue"); + psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_SATURATED, 0, "red"); + psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_PSFSTAR, "yellow"); + psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_SATSTAR, "green"); + 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"); - } + pmVisualAskUser(NULL); return true; } @@ -846,12 +776,6 @@ if (!pmVisualIsVisual()) return true; - if (kapa2 == -1) { - kapa2 = KapaOpenNamedSocket ("kapa", "psphot:psfstars"); - if (kapa2 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } + int myKapa = psphotKapaChannel (3); + if (myKapa == -1) return false; int DX = 64; @@ -898,7 +822,7 @@ 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); + psphotVisualRangeImage (myKapa, psfLogFlux, "psf_mosaic", 0, -2.0, 3.0); + psphotVisualRangeImage (myKapa, funMosaic, "psf_analytical", 1, -10.0, 100.0); + psphotVisualRangeImage (myKapa, resMosaic, "psf_residual", 2, -10.0, 100.0); psFree (psfMosaic); @@ -908,11 +832,5 @@ 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"); - } + pmVisualAskUser(NULL); return true; } @@ -924,12 +842,6 @@ if (!pmVisualIsVisual()) return true; - if (kapa2 == -1) { - kapa2 = KapaOpenNamedSocket ("kapa", "psphot:psfstars"); - if (kapa2 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } + int myKapa = psphotKapaChannel (3); + if (myKapa == -1) return false; // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -1017,10 +929,12 @@ if (Xo == 0) { // place source alone on this row - psphotAddWithTest (source, true, maskVal); // replace source if subtracted + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); - psphotSubWithTest (source, false, maskVal); // remove source (force) + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - psphotSetState (source, false, maskVal); // reset source Add/Sub state to recorded + + if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); Yo += DY; @@ -1032,10 +946,12 @@ Xo = 0; - psphotAddWithTest (source, true, maskVal); // replace source if subtracted + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); - psphotSubWithTest (source, false, maskVal); // remove source (force) + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - psphotSetState (source, false, maskVal); // replace source (has been subtracted) + + if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); Xo = DX; @@ -1044,10 +960,11 @@ } else { // extend this row - psphotAddWithTest (source, true, maskVal); // replace source if subtracted + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); - psphotSubWithTest (source, false, maskVal); // remove source (force) + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - psphotSetState (source, false, maskVal); // replace source (has been subtracted) + if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); Xo += DX; @@ -1056,15 +973,8 @@ } - 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"); - } - + psphotVisualRangeImage (myKapa, outpos, "psfpos", 0, -0.05, 0.95); + psphotVisualRangeImage (myKapa, outsub, "psfsub", 1, -0.05, 0.95); + + pmVisualAskUser(NULL); psFree (outpos); psFree (outsub); @@ -1084,12 +994,6 @@ if (!pmVisualIsVisual()) return true; - if (kapa2 == -1) { - kapa2 = KapaOpenNamedSocket ("kapa", "psphot:images"); - if (kapa2 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } + int myKapa = psphotKapaChannel (3); + if (myKapa == -1) return false; // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -1119,7 +1023,7 @@ pmSource *source = sources->data[i]; - bool keep = false; - keep |= (source->mode & PM_SOURCE_MODE_SATSTAR); - if (!keep) continue; + // only show "real" saturated stars (not defects) + if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;; + if (source->mode & PM_SOURCE_MODE_DEFECT) continue;; // how does this subimage get placed into the output image? @@ -1164,10 +1068,8 @@ pmSource *source = sources->data[i]; - bool keep = false; - if (source->mode & PM_SOURCE_MODE_SATSTAR) { - nSAT ++; - keep = true; - } - if (!keep) continue; + // only show "real" saturated stars (not defects) + if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;; + if (source->mode & PM_SOURCE_MODE_DEFECT) continue;; + nSAT ++; if (Xo + DX > NX) { @@ -1175,7 +1077,8 @@ if (Xo == 0) { // place source alone on this row - psphotAddWithTest (source, true, maskVal); // replace source if subtracted + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); - psphotSetState (source, true, maskVal); // reset source Add/Sub state to recorded + if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); Yo += DY; @@ -1186,7 +1089,9 @@ Yo += dY; Xo = 0; - psphotAddWithTest (source, true, maskVal); // replace source if subtracted + + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); - psphotSetState (source, true, maskVal); // replace source (has been subtracted) + if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); Xo = DX; @@ -1195,7 +1100,8 @@ } else { // extend this row - psphotAddWithTest (source, true, maskVal); // replace source if subtracted + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); - psphotSetState (source, true, maskVal); // replace source (has been subtracted) + if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); Xo += DX; @@ -1204,24 +1110,29 @@ } - 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"); - } - + psphotVisualScaleImage (myKapa, outsat, NULL, "satstar", 2); + + pmVisualAskUser(NULL); psFree (outsat); return true; } +static void plotline (int myKapa, Graphdata *graphdata, float x0, float y0, float x1, float y1) +{ + float x[2], y[2]; + x[0] = x0; + x[1] = x1; + y[0] = y0; + y[1] = y1; + KapaPrepPlot (myKapa, 2, graphdata); + KapaPlotVector (myKapa, 2, x, "x"); + KapaPlotVector (myKapa, 2, y, "y"); +} + 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 + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); int nPts = source->pixels->numRows * source->pixels->numCols; @@ -1240,12 +1151,12 @@ 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] = 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] = 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]); @@ -1255,6 +1166,66 @@ } + // generate model profiles (major and minor axis): + // create a model with theta = 0.0 so major and minor axes are equiv to x and y: + psEllipseShape rawShape, rotShape; + + rawShape.sx = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2; + rawShape.sy = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2; + rawShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY]; + + psEllipseAxes axes = psEllipseShapeToAxes (rawShape, 20.0); + + axes.theta = 0.0; + + rotShape = psEllipseAxesToShape (axes); + + psVector *params = psVectorAlloc(source->modelPSF->params->n, PS_TYPE_F32); + for (int i = 0; i < source->modelPSF->params->n; i++) { + params->data.F32[i] = source->modelPSF->params->data.F32[i]; + } + params->data.F32[PM_PAR_SXX] = rotShape.sx * M_SQRT2; + params->data.F32[PM_PAR_SYY] = rotShape.sy * M_SQRT2; + params->data.F32[PM_PAR_SXY] = rotShape.sxy; + params->data.F32[PM_PAR_XPOS] = 0.0; + params->data.F32[PM_PAR_YPOS] = 0.0; + + psVector *rmod = psVectorAlloc(300, PS_TYPE_F32); + psVector *fmaj = psVectorAlloc(300, PS_TYPE_F32); + psVector *fmin = psVectorAlloc(300, PS_TYPE_F32); + + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + + float r = 0.0; + for (int i = 0; i < rmod->n; i++) { + r = i*0.1; + rmod->data.F32[i] = r; + + coord->data.F32[1] = r; + coord->data.F32[0] = 0.0; + fmaj->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord)); + + coord->data.F32[0] = r; + coord->data.F32[1] = 0.0; + fmin->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord)); + } + psFree (coord); + psFree (params); + + float FWHM_MAJOR = 2.0*source->modelPSF->modelRadius (source->modelPSF->params, 0.5*source->modelPSF->params->data.F32[PM_PAR_I0]); + float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); + if (FWHM_MAJOR < FWHM_MINOR) PS_SWAP (FWHM_MAJOR, FWHM_MINOR); + + psEllipseMoments emoments; + emoments.x2 = source->moments->Mxx; + emoments.xy = source->moments->Mxy; + emoments.y2 = source->moments->Myy; + axes = psEllipseMomentsToAxes (emoments, 20.0); + float MOMENTS_MAJOR = 2.355*axes.major; + float MOMENTS_MINOR = 2.355*axes.minor; + + float logHM = log10(0.5*source->modelPSF->params->data.F32[PM_PAR_I0]); + // reset source Add/Sub state to recorded - psphotSetState (source, state, maskVal); + if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); KapaInitGraph (&graphdata); @@ -1291,4 +1262,33 @@ KapaPlotVector (myKapa, nb, fb->data.F32, "y"); + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 0; + graphdata.size = 0.0; + graphdata.style = 0; + KapaPrepPlot (myKapa, rmod->n, &graphdata); + KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x"); + KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y"); + plotline (myKapa, &graphdata, 0.0, logHM, 30.0, logHM); + plotline (myKapa, &graphdata, 0.5*FWHM_MINOR, 0.0, 0.5*FWHM_MINOR, 5.0); + graphdata.ltype = 1; + plotline (myKapa, &graphdata, 0.5*MOMENTS_MINOR, 0.0, 0.5*MOMENTS_MINOR, 5.0); + graphdata.ltype = 0; + + graphdata.color = KapaColorByName ("green"); + graphdata.ptype = 0; + graphdata.size = 0.0; + graphdata.style = 0; + KapaPrepPlot (myKapa, rmod->n, &graphdata); + KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x"); + KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y"); + plotline (myKapa, &graphdata, 0.5*FWHM_MAJOR, 0.0, 0.5*FWHM_MAJOR, 5.0); + graphdata.ltype = 1; + plotline (myKapa, &graphdata, 0.5*MOMENTS_MAJOR, 0.0, 0.5*MOMENTS_MAJOR, 5.0); + graphdata.ltype = 0; + + for (int i = 0; i < rmod->n; i++) { + rmod->data.F32[i] = log10(rmod->data.F32[i]); + } + // ** loglog ** KapaSelectSection (myKapa, "loglog"); @@ -1299,4 +1299,5 @@ graphdata.ymin = -0.05; graphdata.ymax = +5.05; + graphdata.color = KapaColorByName ("black"); KapaSetLimits (myKapa, &graphdata); @@ -1321,4 +1322,24 @@ KapaPlotVector (myKapa, nb, Rb->data.F32, "x"); KapaPlotVector (myKapa, nb, fb->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 0; + graphdata.size = 0.0; + graphdata.style = 0; + KapaPrepPlot (myKapa, rmod->n, &graphdata); + KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x"); + KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y"); + + graphdata.color = KapaColorByName ("green"); + graphdata.ptype = 0; + graphdata.size = 0.0; + graphdata.style = 0; + KapaPrepPlot (myKapa, rmod->n, &graphdata); + KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x"); + KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y"); + + psFree (rmod); + psFree (fmin); + psFree (fmaj); psFree (rg); @@ -1337,12 +1358,6 @@ if (!pmVisualIsVisual()) return true; - if (kapa3 == -1) { - kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots"); - if (kapa3 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } + int myKapa = psphotKapaChannel (2); + if (myKapa == -1) return false; // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -1351,5 +1366,7 @@ assert (maskVal); - KapaClearPlots (kapa3); + section.bg = KapaColorByName ("none"); // XXX probably should be 'none' + + KapaClearPlots (myKapa); // first section : mag vs CR nSigma section.dx = 1.0; @@ -1359,5 +1376,5 @@ section.name = NULL; psStringAppend (§ion.name, "linlog"); - KapaSetSection (kapa3, §ion); + KapaSetSection (myKapa, §ion); psFree (section.name); @@ -1369,5 +1386,5 @@ section.name = NULL; psStringAppend (§ion.name, "loglog"); - KapaSetSection (kapa3, §ion); + KapaSetSection (myKapa, §ion); psFree (section.name); @@ -1378,5 +1395,5 @@ if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; - psphotVisualPlotRadialProfile (kapa3, source, maskVal); + psphotVisualPlotRadialProfile (myKapa, source, maskVal); // pause and wait for user input: @@ -1388,5 +1405,5 @@ } if (key[0] == 'e') { - KapaClearPlots (kapa3); + KapaClearPlots (myKapa); } if (key[0] == 's') { @@ -1407,10 +1424,11 @@ psEllipseAxes axes; + // XXX skip this for now: it is not very clear + return true; + if (!pmVisualIsVisual()) return true; - if (kapa == -1) { - fprintf (stderr, "kapa not opened, skipping\n"); - return false; - } + int myKapa = psphotKapaChannel (1); + if (myKapa == -1) return false; // note: this uses the Ohana allocation tools: @@ -1485,40 +1503,28 @@ } - KiiLoadOverlay (kapa, overlayE, NoverlayE, "red"); - KiiLoadOverlay (kapa, overlayO, NoverlayO, "yellow"); + KiiLoadOverlay (myKapa, overlayE, NoverlayE, "red"); + KiiLoadOverlay (myKapa, 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; + pmVisualAskUser(NULL); + + return true; +} + +bool psphotVisualShowSourceSize_Single (int myKapa, psArray *sources, pmSourceMode mode, bool keep, float scale, char *color) { + + int Noverlay; KiiOverlay *overlay; - if (!pmVisualIsVisual()) return true; - - if (kapa == -1) { - fprintf (stderr, "kapa not opened, skipping\n"); - return false; - } + psEllipseMoments emoments; + psEllipseAxes axes; // note: this uses the Ohana allocation tools: + ALLOCATE (overlay, KiiOverlay, sources->n); + Noverlay = 0; - NOVERLAY = 100; - ALLOCATE (overlay, KiiOverlay, sources->n); - - // mark CRs with red boxes for (int i = 0; i < sources->n; i++) { @@ -1526,60 +1532,64 @@ 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; + if (mode) { + if (keep) { + if (!(source->mode & mode)) continue; + } else { + if (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 = scale*2.0*axes.major; + overlay[Noverlay].dy = scale*2.0*axes.minor; + overlay[Noverlay].angle = axes.theta * PS_DEG_RAD; 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"); + } + + KiiLoadOverlay (myKapa, overlay, Noverlay, color); 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) { - + return true; +} + +bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources) { + + if (!pmVisualIsVisual()) return true; + + int myKapa = psphotKapaChannel (1); + if (myKapa == -1) return false; + + KiiEraseOverlay (myKapa, "red"); + KiiEraseOverlay (myKapa, "blue"); + KiiEraseOverlay (myKapa, "green"); + KiiEraseOverlay (myKapa, "yellow"); + + psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_EXT_LIMIT | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_SATSTAR, 0, 1.0, "green"); + psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_EXT_LIMIT, 1, 1.0, "blue"); + psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_CR_LIMIT, 1, 1.0, "red"); + psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_DEFECT, 1, 2.0, "red"); + psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_SATSTAR, 1, 1.0, "yellow"); + + fprintf (stdout, "red: CR; blue: EXTENDED; green: PSF-like; yellow: SATSTAR\n"); + pmVisualAskUser(NULL); + return true; +} + +bool psphotVisualPlotSourceSize (psMetadata *recipe, psMetadata *analysis, psArray *sources) { + + bool status; Graphdata graphdata; KapaSection section; @@ -1587,185 +1597,540 @@ if (!pmVisualIsVisual()) return true; - if (kapa3 == -1) { - kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots"); - if (kapa3 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } - - KapaClearPlots (kapa3); + int myKapa = psphotKapaChannel (2); + if (myKapa == -1) return false; + + KapaClearPlots (myKapa); 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 (§ion.name, "a1"); - KapaSetSection (kapa3, §ion); + KapaSetFont (myKapa, "courier", 14); + + section.bg = KapaColorByName ("none"); // XXX probably should be 'none' + + // select the max psfX,Y values for the plot limits + float Xmin = 1000.0, Xmax = 0.0; + float Ymin = 1000.0, Ymax = 0.0; + { + int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS"); + for (int n = 0; n < nRegions; n++) { + + char regionName[64]; + snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n); + psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, 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 X0 = psfX - 10.0*psfdX; + float X1 = psfX + 10.0*psfdX; + float Y0 = psfY - 10.0*psfdY; + float Y1 = psfY + 10.0*psfdY; + + if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); } + if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); } + if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); } + if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); } + } + } + Xmin = PS_MAX(Xmin, -0.1); + Ymin = PS_MAX(Ymin, -0.1); + + // storage vectors for data to be plotted + psVector *xSAT = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *ySAT = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *mSAT = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *sSAT = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + + psVector *xPSF = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *yPSF = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *mPSF = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *sPSF = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + + psVector *xEXT = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *yEXT = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *mEXT = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *sEXT = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + + psVector *xDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *yDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *mDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *sDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + + psVector *xLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *yLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *mLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *sLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + + psVector *xCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *yCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *mCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *sCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + + // construct the vectors + int nSAT = 0; + int nEXT = 0; + int nPSF = 0; + int nDEF = 0; + int nLOW = 0; + int nCR = 0; + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (source->moments == NULL) continue; + + // only plot the measured sources... + if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)) continue; + + if (source->mode & PM_SOURCE_MODE_CR_LIMIT) { + xCR->data.F32[nCR] = source->moments->Mxx; + yCR->data.F32[nCR] = source->moments->Myy; + mCR->data.F32[nCR] = -2.5*log10(source->moments->Sum); + sCR->data.F32[nCR] = source->extNsigma; + nCR++; + } + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + xSAT->data.F32[nSAT] = source->moments->Mxx; + ySAT->data.F32[nSAT] = source->moments->Myy; + mSAT->data.F32[nSAT] = -2.5*log10(source->moments->Sum); + sSAT->data.F32[nSAT] = source->extNsigma; + nSAT++; + } + if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { + xEXT->data.F32[nEXT] = source->moments->Mxx; + yEXT->data.F32[nEXT] = source->moments->Myy; + mEXT->data.F32[nEXT] = -2.5*log10(source->moments->Sum); + sEXT->data.F32[nEXT] = source->extNsigma; + nEXT++; + continue; + } + if (source->mode & PM_SOURCE_MODE_DEFECT) { + xDEF->data.F32[nDEF] = source->moments->Mxx; + yDEF->data.F32[nDEF] = source->moments->Myy; + mDEF->data.F32[nDEF] = -2.5*log10(source->moments->Sum); + sDEF->data.F32[nDEF] = source->extNsigma; + nDEF++; + continue; + } + if (source->errMag > 0.1) { + xLOW->data.F32[nLOW] = source->moments->Mxx; + yLOW->data.F32[nLOW] = source->moments->Myy; + mLOW->data.F32[nLOW] = -2.5*log10(source->moments->Sum); + sLOW->data.F32[nLOW] = source->extNsigma; + nLOW++; + continue; + } + xPSF->data.F32[nPSF] = source->moments->Mxx; + yPSF->data.F32[nPSF] = source->moments->Myy; + mPSF->data.F32[nPSF] = -2.5*log10(source->moments->Sum); + sPSF->data.F32[nPSF] = source->extNsigma; + nPSF++; + } + + xSAT->n = nSAT; + ySAT->n = nSAT; + mSAT->n = nSAT; + sSAT->n = nSAT; + + xPSF->n = nPSF; + yPSF->n = nPSF; + mPSF->n = nPSF; + sPSF->n = nPSF; + + xEXT->n = nEXT; + yEXT->n = nEXT; + mEXT->n = nEXT; + sEXT->n = nEXT; + + xCR->n = nCR; + yCR->n = nCR; + mCR->n = nCR; + sCR->n = nCR; + + xDEF->n = nDEF; + yDEF->n = nDEF; + mDEF->n = nDEF; + sDEF->n = nDEF; + + xLOW->n = nLOW; + yLOW->n = nLOW; + mLOW->n = nLOW; + sLOW->n = nLOW; + + // four sections: MxxMyy, MagMxx, MagMyy, MagSigma + + // first section: MxxMyy + section.dx = 0.75; + section.dy = 0.60; + section.x = 0.00; + section.y = 0.00; + section.name = psStringCopy ("MxxMyy"); + KapaSetSection (myKapa, §ion); psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + graphdata.xmin = Xmin; + graphdata.ymin = Ymin; + graphdata.xmax = Xmax; + graphdata.ymax = Ymax; + KapaSetLimits (myKapa, &graphdata); + + graphdata.padXm = NAN; + graphdata.padYm = NAN; + graphdata.padXp = 0.5; + graphdata.padYp = 0.5; + KapaBox (myKapa, &graphdata); + KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM); + KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nPSF, &graphdata); + KapaPlotVector (myKapa, nPSF, xPSF->data.F32, "x"); + KapaPlotVector (myKapa, nPSF, yPSF->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nEXT, &graphdata); + KapaPlotVector (myKapa, nEXT, xEXT->data.F32, "x"); + KapaPlotVector (myKapa, nEXT, yEXT->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nDEF, &graphdata); + KapaPlotVector (myKapa, nDEF, xDEF->data.F32, "x"); + KapaPlotVector (myKapa, nDEF, yDEF->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 7; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (myKapa, nCR, &graphdata); + KapaPlotVector (myKapa, nCR, xCR->data.F32, "x"); + KapaPlotVector (myKapa, nCR, yCR->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 7; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (myKapa, nSAT, &graphdata); + KapaPlotVector (myKapa, nSAT, xSAT->data.F32, "x"); + KapaPlotVector (myKapa, nSAT, ySAT->data.F32, "y"); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 7; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (myKapa, nLOW, &graphdata); + KapaPlotVector (myKapa, nLOW, xLOW->data.F32, "x"); + KapaPlotVector (myKapa, nLOW, yLOW->data.F32, "y"); + + // second section: MagMyy + section.dx = 0.75; + section.dy = 0.20; + section.x = 0.00; + section.y = 0.80; + section.name = psStringCopy ("MagMyy"); + KapaSetSection (myKapa, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + graphdata.xmin = -17.1; + graphdata.xmax = -6.9; + graphdata.ymin = Ymin; + graphdata.ymax = Ymax; + KapaSetLimits (myKapa, &graphdata); + + graphdata.padXm = 0.5; + graphdata.padYm = NAN; + graphdata.padXp = NAN; + graphdata.padYp = 0.5; + strcpy (graphdata.labels, "0210"); + KapaBox (myKapa, &graphdata); + KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_XP); + KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nPSF, &graphdata); + KapaPlotVector (myKapa, nPSF, mPSF->data.F32, "x"); + KapaPlotVector (myKapa, nPSF, yPSF->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nEXT, &graphdata); + KapaPlotVector (myKapa, nEXT, mEXT->data.F32, "x"); + KapaPlotVector (myKapa, nEXT, yEXT->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nDEF, &graphdata); + KapaPlotVector (myKapa, nDEF, mDEF->data.F32, "x"); + KapaPlotVector (myKapa, nDEF, yDEF->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 7; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (myKapa, nCR, &graphdata); + KapaPlotVector (myKapa, nCR, mCR->data.F32, "x"); + KapaPlotVector (myKapa, nCR, yCR->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 7; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (myKapa, nSAT, &graphdata); + KapaPlotVector (myKapa, nSAT, mSAT->data.F32, "x"); + KapaPlotVector (myKapa, nSAT, ySAT->data.F32, "y"); + + // third section: MagMxx + section.dx = 0.25; + section.dy = 0.60; + section.x = 0.75; + section.y = 0.00; + section.name = psStringCopy ("MagMxx"); + KapaSetSection (myKapa, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + graphdata.xmin = Xmin; + graphdata.xmax = Xmax; + graphdata.ymin = -6.9; + graphdata.ymax = -17.1; + KapaSetLimits (myKapa, &graphdata); + + graphdata.padXm = NAN; + graphdata.padYm = 0.5; + graphdata.padXp = 0.5; + graphdata.padYp = NAN; + strcpy (graphdata.labels, "2001"); + KapaBox (myKapa, &graphdata); + KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM); + KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_YP); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nPSF, &graphdata); + KapaPlotVector (myKapa, nPSF, xPSF->data.F32, "x"); + KapaPlotVector (myKapa, nPSF, mPSF->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nEXT, &graphdata); + KapaPlotVector (myKapa, nEXT, xEXT->data.F32, "x"); + KapaPlotVector (myKapa, nEXT, mEXT->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nDEF, &graphdata); + KapaPlotVector (myKapa, nDEF, xDEF->data.F32, "x"); + KapaPlotVector (myKapa, nDEF, mDEF->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 7; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (myKapa, nCR, &graphdata); + KapaPlotVector (myKapa, nCR, xCR->data.F32, "x"); + KapaPlotVector (myKapa, nCR, mCR->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 7; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (myKapa, nSAT, &graphdata); + KapaPlotVector (myKapa, nSAT, xSAT->data.F32, "x"); + KapaPlotVector (myKapa, nSAT, mSAT->data.F32, "y"); + + // fourth section: MagSigma + section.dx = 0.75; + section.dy = 0.20; + section.x = 0.00; + section.y = 0.60; + section.name = psStringCopy ("MagSigma"); + KapaSetSection (myKapa, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + graphdata.xmax = -6.9; + graphdata.xmin = -17.1; + graphdata.ymin = -20.1; + graphdata.ymax = +20.1; + KapaSetLimits (myKapa, &graphdata); + + graphdata.padXm = 0.5; + graphdata.padYm = NAN; + graphdata.padXp = 0.5; + graphdata.padYp = 0.5; + strcpy (graphdata.labels, "0100"); + KapaBox (myKapa, &graphdata); + // KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_XM); + KapaSendLabel (myKapa, "EXT&ss&c", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nPSF, &graphdata); + KapaPlotVector (myKapa, nPSF, mPSF->data.F32, "x"); + KapaPlotVector (myKapa, nPSF, sPSF->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nEXT, &graphdata); + KapaPlotVector (myKapa, nEXT, mEXT->data.F32, "x"); + KapaPlotVector (myKapa, nEXT, sEXT->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nDEF, &graphdata); + KapaPlotVector (myKapa, nDEF, mDEF->data.F32, "x"); + KapaPlotVector (myKapa, nDEF, sDEF->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 7; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (myKapa, nCR, &graphdata); + KapaPlotVector (myKapa, nCR, mCR->data.F32, "x"); + KapaPlotVector (myKapa, nCR, sCR->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 7; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (myKapa, nSAT, &graphdata); + KapaPlotVector (myKapa, nSAT, mSAT->data.F32, "x"); + KapaPlotVector (myKapa, nSAT, sSAT->data.F32, "y"); + + // draw N circles to outline the clumps + { + KapaSelectSection (myKapa, "MxxMyy"); + + // 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, analysis, "PSF.CLUMP.NREGIONS"); + float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.style = 0; + + graphdata.xmin = Xmin; + graphdata.ymin = Ymin; + graphdata.xmax = Xmax; + graphdata.ymax = Ymax; + KapaSetLimits (myKapa, &graphdata); + + for (int n = 0; n < nRegions; n++) { + + char regionName[64]; + snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n); + psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, 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 (myKapa, xLimit->n, &graphdata); + KapaPlotVector (myKapa, xLimit->n, xLimit->data.F32, "x"); + KapaPlotVector (myKapa, yLimit->n, yLimit->data.F32, "y"); + } + psFree (xLimit); + psFree (yLimit); + } + + psFree (xSAT); + psFree (ySAT); + psFree (mSAT); + psFree (sSAT); + + psFree (xEXT); + psFree (yEXT); + psFree (mEXT); + psFree (sEXT); + + psFree (xPSF); + psFree (yPSF); + psFree (mPSF); + psFree (sPSF); + + psFree (xDEF); + psFree (yDEF); + psFree (mDEF); + psFree (sDEF); + + psFree (xLOW); + psFree (yLOW); + psFree (mLOW); + psFree (sLOW); + + psFree (xCR); + psFree (yCR); + psFree (mCR); + psFree (sCR); + + pmVisualAskUser(NULL); + return true; +} + +bool psphotVisualShowResidualImage (pmReadout *readout) { + + if (!pmVisualIsVisual()) return true; + + int myKapa = psphotKapaChannel (1); + if (myKapa == -1) return false; + + psphotVisualScaleImage (myKapa, readout->image, readout->mask, "resid", 1); + + pmVisualAskUser(NULL); + return true; +} + +bool psphotVisualPlotApResid (psArray *sources, float mean, float error) { + + Graphdata graphdata; + float lineX[2], lineY[2]; + + if (!pmVisualIsVisual()) return true; + + int myKapa = psphotKapaChannel (2); + if (myKapa == -1) return false; + + KapaClearPlots (myKapa); + 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->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 (§ion.name, "a2"); - KapaSetSection (kapa3, §ion); - 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 (!pmVisualIsVisual()) return true; - - if (kapa == -1) { - kapa = KapaOpenNamedSocket ("kapa", "psphot:images"); - if (kapa == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(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 (!pmVisualIsVisual()) return true; - - if (kapa3 == -1) { - kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots"); - if (kapa3 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } - - KapaClearPlots (kapa3); - KapaInitGraph (&graphdata); - - psVector *x = psVectorAllocEmpty (sources->n, PS_TYPE_F32); - psVector *y = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *dy = psVectorAllocEmpty (sources->n, PS_TYPE_F32); graphdata.xmin = +32.0; @@ -1785,4 +2150,5 @@ x->data.F32[n] = source->psfMag; y->data.F32[n] = source->apMag - source->psfMag; + dy->data.F32[n] = source->errMag; graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]); graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]); @@ -1792,5 +2158,5 @@ n++; } - x->n = y->n = n; + x->n = y->n = dy->n = n; float range; @@ -1802,11 +2168,16 @@ 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); + // XXX test + graphdata.xmin = -17.0; + graphdata.xmax = -9.0; + graphdata.ymin = -0.31; + graphdata.ymax = +0.31; + + KapaSetLimits (myKapa, &graphdata); + + KapaSetFont (myKapa, "helvetica", 14); + KapaBox (myKapa, &graphdata); + KapaSendLabel (myKapa, "PSF Mag", KAPA_LABEL_XM); + KapaSendLabel (myKapa, "Ap Mag - PSF Mag", KAPA_LABEL_YM); graphdata.color = KapaColorByName ("black"); @@ -1814,18 +2185,187 @@ 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"); + graphdata.etype |= 0x01; + KapaPrepPlot (myKapa, n, &graphdata); + KapaPlotVector (myKapa, n, x->data.F32, "x"); + KapaPlotVector (myKapa, n, y->data.F32, "y"); + KapaPlotVector (myKapa, n, dy->data.F32, "dym"); + KapaPlotVector (myKapa, n, dy->data.F32, "dyp"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 0; + graphdata.etype = 0; + lineX[0] = graphdata.xmin; + lineX[1] = graphdata.xmax; + lineY[0] = lineY[1] = mean; + KapaPrepPlot (myKapa, 2, &graphdata); + KapaPlotVector (myKapa, 2, lineX, "x"); + KapaPlotVector (myKapa, 2, lineY, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 0; + graphdata.etype = 0; + lineX[0] = graphdata.xmin; + lineX[1] = graphdata.xmax; + lineY[0] = lineY[1] = mean + error; + KapaPrepPlot (myKapa, 2, &graphdata); + KapaPlotVector (myKapa, 2, lineX, "x"); + KapaPlotVector (myKapa, 2, lineY, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 0; + graphdata.etype = 0; + lineX[0] = graphdata.xmin; + lineX[1] = graphdata.xmax; + lineY[0] = lineY[1] = mean - error; + KapaPrepPlot (myKapa, 2, &graphdata); + KapaPlotVector (myKapa, 2, lineX, "x"); + KapaPlotVector (myKapa, 2, lineY, "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"); - } + psFree (dy); + + pmVisualAskUser(NULL); + return true; +} + +bool psphotVisualPlotChisq (psArray *sources) { + + Graphdata graphdata; + + if (!pmVisualIsVisual()) return true; + + int myKapa = psphotKapaChannel (2); + if (myKapa == -1) return false; + + KapaClearPlots (myKapa); + 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; + + FILE *f = fopen ("chisq.dat", "w"); + + // 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 (!source->moments) continue; + if (!isfinite(source->moments->Sum)) continue; + if (!source->modelPSF) continue; + if (!isfinite(source->modelPSF->chisq)) continue; + + x->data.F32[n] = -2.5*log10(source->moments->Sum); + y->data.F32[n] = source->modelPSF->chisq / source->modelPSF->nDOF; + 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]); + + fprintf (f, "%d %d %f %f\n", i, n, x->data.F32[n], y->data.F32[n]); + + n++; + } + x->n = y->n = n; + fclose (f); + + 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 test + graphdata.xmin = -17.0; + graphdata.xmax = -3.0; + graphdata.ymin = -0.1; + graphdata.ymax = +10.1; + + KapaSetLimits (myKapa, &graphdata); + + KapaSetFont (myKapa, "helvetica", 14); + KapaBox (myKapa, &graphdata); + KapaSendLabel (myKapa, "PSF Mag", KAPA_LABEL_XM); + KapaSendLabel (myKapa, "ChiSq", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 2; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, n, &graphdata); + KapaPlotVector (myKapa, n, x->data.F32, "x"); + KapaPlotVector (myKapa, n, y->data.F32, "y"); + + psFree (x); + psFree (y); + + pmVisualAskUser(NULL); + return true; +} + +bool psphotVisualShowPetrosians (psArray *sources) { + + int Noverlay, NOVERLAY; + KiiOverlay *overlay; + + if (!pmVisualIsVisual()) return true; + + int kapa = psphotKapaChannel (1); + if (kapa == -1) return false; + + Noverlay = 0; + NOVERLAY = 100; + ALLOCATE (overlay, KiiOverlay, NOVERLAY); + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + if (!source) continue; + if (!source->extpars) continue; + if (!source->extpars->petProfile) continue; + + float petrosianRadius = source->extpars->petrosianRadius; + psEllipseAxes *axes = &source->extpars->axes; + + overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + overlay[Noverlay].x = source->peak->xf; + overlay[Noverlay].y = source->peak->yf; + overlay[Noverlay].dx = 1.0*petrosianRadius; + overlay[Noverlay].dy = 1.0*petrosianRadius*axes->minor/axes->major; + overlay[Noverlay].angle = axes->theta * PS_DEG_RAD; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + + // overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + // overlay[Noverlay].x = source->peak->xf; + // overlay[Noverlay].y = source->peak->yf; + // overlay[Noverlay].dx = 2.0*petrosianRadius; + // overlay[Noverlay].dy = 2.0*petrosianRadius*axes->minor/axes->major; + // overlay[Noverlay].angle = axes->theta * PS_DEG_RAD; + // overlay[Noverlay].text = NULL; + // Noverlay ++; + // CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + } + + KiiLoadOverlay (kapa, overlay, Noverlay, "red"); + FREE (overlay); + + pmVisualAskUser(NULL); return true; } @@ -1850,2 +2390,40 @@ # endif + +# 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 (myKapa, dX->n, &graphdata); +KapaPlotVector (myKapa, dX->n, dX->data.F32, "x"); +KapaPlotVector (myKapa, vX->n, vX->data.F32, "y"); + +psFree (nX); +psFree (dX); +psFree (vX); + +psFree (nY); +psFree (dY); +psFree (vY); + +# endif +