Index: trunk/psModules/src/astrom/pmAstrometryObjects.c =================================================================== --- trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 21387) +++ trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 21422) @@ -8,6 +8,6 @@ * @author EAM, IfA * -* @version $Revision: 1.44 $ $Name: not supported by cvs2svn $ -* @date $Date: 2009-02-01 21:45:33 $ +* @version $Revision: 1.45 $ $Name: not supported by cvs2svn $ +* @date $Date: 2009-02-09 21:25:20 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -32,4 +32,6 @@ #include "pmHDU.h" #include "pmFPA.h" +#include "pmFPAExtent.h" +#include "pmFPAfile.h" #include "pmAstrometryObjects.h" #include "pmKapaPlots.h" @@ -658,5 +660,5 @@ } - // XXX this function is crashing + // XXX this function is crashing // pmAstromVisualPlotGridMatch(raw, ref, gridNP, stats->offset.x, stats->offset.y, maxOffpix, Scale, Offset); Index: trunk/psModules/src/astrom/pmAstrometryVisual.c =================================================================== --- trunk/psModules/src/astrom/pmAstrometryVisual.c (revision 21387) +++ trunk/psModules/src/astrom/pmAstrometryVisual.c (revision 21422) @@ -1,8 +1,2 @@ -/** Diagnostic plots called from within pmAstrometry routines. - * @author Chris Beaumont - * @date October 2008 - */ - -/* Include Files */ #ifdef HAVE_CONFIG_H #include @@ -16,16 +10,15 @@ #include -#include "pmKapaPlots.h" - #include "pmHDU.h" #include "pmFPA.h" +#include "pmFPAfile.h" #include "pmAstrometryObjects.h" +#include "pmAstrometryVisual.h" +#include "pmFPAExtent.h" # if (HAVE_KAPA) # include - -# define KAPAX 700 -# define KAPAY 700 - +# include "pmKapaPlots.h" +#include "pmVisual.h" //variables to determine when things are plotted @@ -33,17 +26,26 @@ static bool plotGridMatch = true; static bool plotTweak = true; +static bool plotRawStars = true; +static bool plotRefStars = false; +static bool plotLumFunc = true; +static bool plotRemoveClumps = true; +static bool plotOneChipFit = true; +static bool plotFixChips = true; +static bool plotAstromGuessCheck = true; +static bool plotMosaicMatches = true; +static bool plotCommonScale = true; +static bool plotMosaicOneChip = true; // variables to store plotting window indices static int kapa = -1; - - -/* Utility Routine Prototypes */ -bool pmAstromVisualScaleGraphdata(Graphdata *graphdata, psVector *xVec, psVector *yVec, bool clip); +static int kapa2 = -1; + +//helper prototypes +bool residPlot (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe, + char *title); /* Initialization Routines */ - -/** start or stop plotting */ bool pmAstromSetVisual (bool mode) { isVisual = mode; @@ -52,68 +54,762 @@ -/** open name, and resize a window if necessary */ -bool pmAstromVisualInitWindow (int *kapid, char *name) { - if (*kapid == -1) { - *kapid = KapaOpenNamedSocket("kapa", name); - if (*kapid == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled for pmAstrom\n"); - isVisual = false; - return false; - } - KapaResize (*kapid, KAPAX, KAPAY); - } - return true; -} - - -/** ask the user how to proceed */ -bool pmAstromVisualAskUser(bool *plotflag) -{ - char key[10]; - fprintf (stdout, "[c]ontinue? [s]kip the rest of these plots? [a]bort all visual plots?"); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } - if (key[0] == 's') { - *plotflag = false; - } - if (key[0] == 'a') { - pmAstromSetVisual(false); - } - return true; -} - - -/** destroy windows at the end of a run*/ bool pmAstromVisualClose() { if(kapa != -1) KiiClose(kapa); + if(kapa2 != -1) + KiiClose(kapa2); return true; } -/* plotting routines */ - -/** - * Plot the offset between every pair of reference and raw source locations. The peak of this - * distribution nominally gives the offset, scale difference, and rotation of the two catalogs. - * Overplots the location of this peak as determined by pmAstromGridAngle, as well as some profiles - * along horizontal and vertical cuts through this peak. - */ -bool pmAstromVisualPlotGridMatch (const psArray *raw, ///< raw stars - const psArray *ref, ///< reference stars - psImage *gridNP, ///< a 2D histogram of raw-ref star distances - double offsetX, ///< The X location (FP coordinates) of the peak of gridNP - double offsetY, ///< the Y location (FP coordinates) of the peak of gridNP - double maxOffpix, ///< The half-width of gridNP in FP coordinates - double Scale, ///< The pixel size of gridNP in histogram-bin-coordinates - double Offset ///< The (x,y) location (histogram-bin coordinates) of the FP point (0,0) in gridNP +/* Plotting Routines */ +bool pmAstromVisualPlotRawStars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe) +{ + // make sure we want to plot this + if (!plotRawStars || !isVisual) return true; + + //set up plot region + if (!pmVisualInitWindow (&kapa, "psastro:plots")){ + isVisual = false; + return false; + } + Graphdata graphdata; + KapaSection section; + + KapaInitGraph (&graphdata); + KapaClearPlots (kapa); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 7; + graphdata.size = 0.5; + graphdata.style = 2; + + section.dx = 0.4; + section.dy = 0.4; + + //initialize and populate plotting vectors + bool status = false; + float iMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MIN"); + float iMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MAX"); + + psVector *xVec = psVectorAlloc (rawstars->n, PS_TYPE_F32); + psVector *yVec = psVectorAlloc (rawstars->n, PS_TYPE_F32); + psVector *zVec = psVectorAlloc (rawstars->n, PS_TYPE_F32); + + section.x = 0.0; + section.y = 0.5; + section.name = NULL; + psStringAppend (§ion.name, "a0"); + KapaSetSection (kapa, §ion); + psFree (section.name); + + //Chip coordinates + int n = 0; + for (int i = 0; i < rawstars->n; i++) { + pmAstromObj *raw = rawstars->data[i]; + if (!isfinite(raw->Mag)) continue; + if (raw->Mag < iMagMin) continue; + if (raw->Mag > iMagMax) continue; + + xVec->data.F32[n] = raw->chip->x; + yVec->data.F32[n] = raw->chip->y; + zVec->data.F32[n] = raw->Mag; + n++; + } + xVec->n = yVec->n = zVec->n = n; + + KapaSendLabel(kapa, "Chip", KAPA_LABEL_XP); + KapaSendLabel(kapa, "X", KAPA_LABEL_XM); + KapaSendLabel(kapa, "Y", KAPA_LABEL_YM); + pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false); + + //Focal Plane Coordinates + section.x = 0.5; + section.y = 0.5; + section.name = NULL; + psStringAppend (§ion.name, "a1"); + KapaSetSection (kapa, §ion); + psFree (section.name); + + n = 0; + for (int i = 0; i < rawstars->n; i++) { + pmAstromObj *raw = rawstars->data[i]; + if (!isfinite(raw->Mag)) continue; + if (raw->Mag < iMagMin) continue; + if (raw->Mag > iMagMax) continue; + + xVec->data.F32[n] = raw->FP->x; + yVec->data.F32[n] = raw->FP->y; + zVec->data.F32[n] = raw->Mag; + n++; + } + xVec->n = yVec->n = zVec->n = n; + + KapaSendLabel (kapa, "Focal Plane", KAPA_LABEL_XP); + KapaSendLabel (kapa, "L", KAPA_LABEL_XM); + KapaSendLabel (kapa, "M", KAPA_LABEL_YM); + pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false); + + // Tangent Plane Coordinates + section.x = 0.0; + section.y = 0.0; + section.name = NULL; + psStringAppend (§ion.name, "a2"); + KapaSetSection (kapa, §ion); + psFree (section.name); + + n = 0; + for (int i = 0; i < rawstars->n; i++) { + pmAstromObj *raw = rawstars->data[i]; + if (!isfinite(raw->Mag)) continue; + if (raw->Mag < iMagMin) continue; + if (raw->Mag > iMagMax) continue; + + xVec->data.F32[n] = raw->TP->x; + yVec->data.F32[n] = raw->TP->y; + zVec->data.F32[n] = raw->Mag; + n++; + } + xVec->n = yVec->n = zVec->n = n; + + KapaSendLabel (kapa, "Tangential Plane", KAPA_LABEL_XP); + KapaSendLabel (kapa, "P", KAPA_LABEL_XM); + KapaSendLabel (kapa, "Q", KAPA_LABEL_YM); + pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false); + + //sky coordinates + section.x = 0.5; + section.y = 0.0; + section.name = NULL; + psStringAppend (§ion.name, "a3"); + KapaSetSection (kapa, §ion); + psFree (section.name); + + n = 0; + for (int i = 0; i < rawstars->n; i++) { + pmAstromObj *raw = rawstars->data[i]; + if (!isfinite(raw->Mag)) continue; + if (raw->Mag < iMagMin) continue; + if (raw->Mag > iMagMax) continue; + + xVec->data.F32[n] = DEG_RAD*raw->sky->r; + yVec->data.F32[n] = DEG_RAD*raw->sky->d; + zVec->data.F32[n] = raw->Mag; + n++; + } + xVec->n = yVec->n = zVec->n = n; + + KapaSendLabel (kapa, "Sky", KAPA_LABEL_XP); + KapaSendLabel(kapa, "RA", KAPA_LABEL_XM); + KapaSendLabel(kapa, "Dec", KAPA_LABEL_YM); + pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false); + + // flip x (East increase to left) + SWAP (graphdata.xmin, graphdata.xmax); + KapaSetLimits (kapa, &graphdata); + + // plot label + section.x = 0.0; + section.y = 0.0; + section.dx = .95; + section.dy = .95; + section.name = NULL; + psStringAppend (§ion.name, "a5"); + KapaSetSection (kapa, §ion); + KapaSendLabel (kapa, "Raw Star Selection - Initial Astrometry", KAPA_LABEL_XP); + psFree (section.name); + + // pause and wait for user input: + pmVisualAskUser(&plotRawStars, &isVisual); + + psFree (xVec); + psFree (yVec); + psFree (zVec); + return true; +} + + +bool pmAstromVisualPlotRefStars (psArray *refstars, psMetadata *recipe) +{ + //make sure we want to plot this + if (!isVisual || !plotRefStars) return true; + + //set up plotting variables + if (!pmVisualInitWindow (&kapa, "psastro:plots")) { + isVisual = false; + return false; + } + + Graphdata graphdata; + KapaInitGraph (&graphdata); + KapaClearSections (kapa); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 7; + graphdata.size = 0.5; + graphdata.style = 2; + + //initialize and populate plot vectors + bool status = false; + float rMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MIN"); + float rMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MAX"); + + psVector *xVec = psVectorAlloc (refstars->n, PS_TYPE_F32); + psVector *yVec = psVectorAlloc (refstars->n, PS_TYPE_F32); + psVector *zVec = psVectorAlloc (refstars->n, PS_TYPE_F32); + + int n = 0; + for (int i = 0; i < refstars->n; i++) { + pmAstromObj *ref = refstars->data[i]; + if (!isfinite(ref->Mag)) continue; + if (ref->Mag > rMagMax) continue; + if (ref->Mag < rMagMin) continue; + + xVec->data.F32[n] = DEG_RAD*ref->sky->r; + yVec->data.F32[n] = DEG_RAD*ref->sky->d; + zVec->data.F32[n] = ref->Mag; + n++; + } + xVec->n = yVec->n = zVec->n = n; + + KapaSendLabel (kapa, "RA", KAPA_LABEL_XM); + KapaSendLabel (kapa, "Dec", KAPA_LABEL_YM); + KapaSendLabel (kapa, "Reference Stars", KAPA_LABEL_XP); + pmVisualTriplePlot(kapa, &graphdata, xVec, yVec, zVec, false); + + // flip x (East increase to left) + SWAP (graphdata.xmin, graphdata.xmax); + KapaSetLimits (kapa, &graphdata); + + // pause and wait for user input: + pmVisualAskUser(&plotRefStars, &isVisual); + + psFree (xVec); + psFree (yVec); + psFree (zVec); + return true; +} + + +bool pmAstromVisualPlotLuminosityFunction (psVector *lnMag, // Log(n) for each magnitude bin + psVector *Mag, // magnitude bins + pmLumFunc *lumFunc,// Fit to the reference star luminosity function + pmLumFunc *rawFunc // Fit to the raw star luminoisty function + ) +{ + + // make sure we want to plot this + if ( !isVisual || !plotLumFunc ) return true; + + //set up plotting variables + if ( !pmVisualInitWindow (&kapa, "psastro:plots")){ + isVisual = false; + return false; + } + + Graphdata graphdata; + KapaSection section1 = {"s1", .1, .1, .85, .35}; + KapaSection section2 = {"s2", .1, .5, .85, .35}; + KapaSection section; + if(rawFunc == NULL) + section = section1; + else + section = section2; + if (rawFunc == NULL) + KapaClearPlots(kapa); + KapaInitGraph(&graphdata); + + //Determine Plot Limits + pmVisualScaleGraphdata(&graphdata, Mag, lnMag, false); + + //Make a line for the fit + float x[2] = {graphdata.xmin, graphdata.xmax}; + float y[2] = {lumFunc->offset + x[0] * lumFunc->slope, + lumFunc->offset + x[1] * lumFunc->slope}; + + //Plot Data + KapaSetSection(kapa, §ion); + KapaSetLimits(kapa, &graphdata); + + KapaSetFont (kapa, "helvetica", 14); + KapaBox (kapa, &graphdata); + KapaSendLabel (kapa, "Magnitude", KAPA_LABEL_XM); + KapaSendLabel (kapa, "Log(N)", KAPA_LABEL_YM); + if (rawFunc == NULL) + KapaSendLabel (kapa, "Raw Star Luminosity Function", KAPA_LABEL_XP); + else + KapaSendLabel (kapa, + "Reference Star Luminosity Function, Shifted Raw Fit, and Cutoff", + KAPA_LABEL_XP); + graphdata.color=KapaColorByName("black"); + graphdata.style = 1; + KapaPrepPlot (kapa, lnMag->n, &graphdata); + KapaPlotVector(kapa, lnMag->n, Mag->data.F32, "x"); + KapaPlotVector(kapa, lnMag->n, lnMag->data.F32, "y"); + + //Overplot fit + graphdata.style=0; + KapaPrepPlot(kapa,2,&graphdata); + KapaPlotVector(kapa, 2, x, "x"); + KapaPlotVector(kapa, 2, y, "y"); + + //If rawFunc was supplied, overplot the raw star fit + cutoff + if( rawFunc != NULL) { + double mRef = 0.5*(lumFunc->mMin + lumFunc->mMax); + double logRho = mRef * lumFunc->slope + lumFunc->offset; + double mRaw = (logRho - rawFunc->offset) / rawFunc->slope; + double deltaM = mRef - mRaw; + double mRefMax = rawFunc->mMax + deltaM; + + float xraw[2] = {rawFunc->mMin + deltaM, rawFunc->mMax + deltaM}; + float yraw[2] = {rawFunc->offset + (rawFunc->slope) * rawFunc->mMin, + rawFunc->offset + (rawFunc->slope) * rawFunc->mMax}; + float x[2] = {mRefMax, mRefMax}; + float y[2] = {graphdata.ymin, graphdata.ymax}; + graphdata.color= KapaColorByName("red"); + KapaPrepPlot(kapa, 2, &graphdata); + KapaPlotVector(kapa, 2, x, "x"); + KapaPlotVector(kapa, 2, y, "y"); + KapaPrepPlot (kapa, 2, &graphdata); + KapaPlotVector (kapa, 2, xraw, "x"); + KapaPlotVector (kapa, 2, yraw, "y"); + + // pause and wait for user input: + pmVisualAskUser (&plotLumFunc, &isVisual); + } + return true; +} // end of pmAstromVisualPlotLuminosityFunction + + +bool pmAstromVisualPlotRemoveClumps (psArray *input, // Array containing the field stars + psImage *count, // A 2D histogram of the field star distribution + int scale, // The pixel size of the histogram + float limit // The minimum numuber of stars in a bin flagged as a clump + ) +{ + + //make sure we want to plot this + if (!isVisual || !plotRemoveClumps) return true; + + //set up plot variables + if ( !pmVisualInitWindow (&kapa, "psastro:plots")) { + isVisual = false; + return false; + } + + KapaSection section; + Graphdata graphdata; + section.x = 0; + section.dx = .9; + section.y = 0; + section.dy = .9; + section.name = NULL; + psStringAppend( §ion.name, "a0"); + KapaInitGraph(&graphdata); + KapaSetSection(kapa, §ion); + psFree(section.name); + + graphdata.ptype = 7; + graphdata.size = 0.5; + graphdata.style = 2; + graphdata.color = KapaColorByName ("black"); + KapaClearPlots(kapa); + + //set up plot vectors + float Xmin = +FLT_MAX; + float Xmax = -FLT_MAX; + float Ymin = +FLT_MAX; + float Ymax = -FLT_MAX; + psVector *xVec = psVectorAlloc (input->n, PS_TYPE_F32); + psVector *yVec = psVectorAlloc (input->n, PS_TYPE_F32); + + //determine boundaries for histogram bin calculation + int n = 0; + for (int i=0; i< input->n; i++) { + pmAstromObj *obj = (pmAstromObj *)input->data[i]; + if (!isfinite(obj->FP->x)) continue; + if (!isfinite(obj->FP->y)) continue; + xVec->data.F32[n] = obj->FP->x; + yVec->data.F32[n] = obj->FP->y; + Xmin = PS_MIN (Xmin, xVec->data.F32[n]); + Xmax = PS_MAX (Xmax, xVec->data.F32[n]); + Ymin = PS_MIN (Ymin, yVec->data.F32[n]); + Ymax = PS_MAX (Ymax, yVec->data.F32[n]); + n++; + } + xVec->n = yVec->n = n; + + //plot stars + graphdata.xmax = Xmax; + graphdata.xmin = Xmin; + graphdata.ymax = Ymax; + graphdata.ymin = Ymin; + KapaSetLimits (kapa, &graphdata); + KapaSetFont (kapa, "helvetica", 14); + + KapaBox (kapa, &graphdata); + + KapaSendLabel (kapa, "L (pixels)", KAPA_LABEL_XM); + KapaSendLabel (kapa, "M (pixels)", KAPA_LABEL_YM); + KapaSendLabel (kapa, "Regions Flagged as Clumps (Red Boxes)", + KAPA_LABEL_XP); + + KapaPrepPlot (kapa, xVec->n, &graphdata); + KapaPlotVector (kapa, xVec->n, xVec->data.F32, "x"); + KapaPlotVector (kapa, yVec->n, yVec->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.style = 1; + + //overplot clumpy regions excluded from analysis + for(int i = 0; i < count->numCols; i++) { + for (int j = 0; j < count->numRows; j++) { + if(count->data.U32[j][i] <= limit) continue; //not a clump + float Xbot = (i - 5) * scale + Xmin; + float Ybot = (j - 5) * scale + Ymin; + if(Xbot < graphdata.xmin || Xbot > graphdata.xmax || + Ybot < graphdata.ymin || Ybot > graphdata.ymax) continue; + float x[5] = {Xbot, Xbot + scale, Xbot + scale, Xbot, Xbot}; + float y[5] = {Ybot, Ybot, Ybot + scale, Ybot + scale, Ybot}; + KapaPrepPlot (kapa, 5, &graphdata); + KapaPlotVector (kapa, 5, x, "x"); + KapaPlotVector (kapa, 5, y, "y"); + } + } + + //ask for user input and finish + pmVisualAskUser (&plotRemoveClumps, &isVisual); + psFree (xVec); + psFree (yVec); + + return true; +} + + +bool pmAstromVisualPlotOneChipFit (psArray *rawstars, // stars detected in the image + psArray *refstars, // reference stars over the same region + psArray *match, // contains which rawstars match to which refstars + psMetadata *recipe // data reduction recipe ) +{ + + //make sure we want to plot this + if (!isVisual || !plotOneChipFit) + return true; + + if(!pmVisualInitWindow(&kapa, "psastro:plots") || !pmVisualInitWindow(&kapa2, "psastro:plots")) { + isVisual = false; + return false; + } + + //plot the residuals + if (!residPlot(rawstars, refstars, match, recipe, "Single Chip Fit Residuals (Chip Coordinates)")) { + plotOneChipFit = false; + return false; + } + + //ask for user input and finish + pmVisualAskUser(&plotOneChipFit, &isVisual); + return true; +} + + +bool pmAstromVisualPlotFixChips (pmFPAfile *input, // focal plane array file + psVector *xOld, // old X location of chip cornerss + psVector *yOld // old Y location of chip corners + ) +{ + //make sure we want to plot this + if(!isVisual || !plotFixChips) return true; + + if(!pmVisualInitWindow(&kapa, "psastro:plots")) { + isVisual = false; + return false; + } + + KapaSection section = {"s1", .05, .05, .9, .9}; + Graphdata graphdata; + KapaInitGraph( &graphdata); + KapaClearPlots (kapa); + graphdata.ptype = 2; + graphdata.style = 2; + + psVector *xNew = psVectorAllocEmpty (xOld->n, PS_TYPE_F32); + psVector *yNew = psVectorAllocEmpty (yOld->n, PS_TYPE_F32); + + // copy of the code in psastroFixChips that generated xOld, yOld, but for xNew, yNew + pmFPAview *view = pmFPAviewAlloc (0); + + pmChip *obsChip = NULL; + while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { + if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; } + + psRegion *region = pmChipPixels(obsChip); + psPlane ptCP, ptFP; + + ptCP.x = region->x0; ptCP.y = region->y0; + psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); + psVectorAppend (xNew, ptFP.x); + psVectorAppend (yNew, ptFP.y); + + ptCP.x = region->x0; ptCP.y = region->y1; + psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); + psVectorAppend (xNew, ptFP.x); + psVectorAppend (yNew, ptFP.y); + + ptCP.x = region->x1; ptCP.y = region->y1; + psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); + psVectorAppend (xNew, ptFP.x); + psVectorAppend (yNew, ptFP.y); + + ptCP.x = region->x1; ptCP.y = region->y0; + psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); + psVectorAppend (xNew, ptFP.x); + psVectorAppend (yNew, ptFP.y); + + psFree (region); + } + + //set up graph + pmVisualScaleGraphdata(&graphdata, xOld, yOld, true); + pmVisualInitGraph(kapa, §ion, &graphdata); + KapaSendLabel (kapa, "L (FP)", KAPA_LABEL_XM); + KapaSendLabel (kapa, "M (FP)", KAPA_LABEL_YM); + KapaSendLabel (kapa, "Chip corners before (black) and after (red) FixChips", KAPA_LABEL_XP); + KapaPrepPlot (kapa, xOld->n, &graphdata); + KapaPlotVector (kapa, xOld->n, xOld->data.F32, "x"); + KapaPlotVector (kapa, xOld->n, yOld->data.F32, "y"); + + graphdata.ptype = 1; + graphdata.color = KapaColorByName("red"); + KapaPrepPlot (kapa, xNew->n, &graphdata); + KapaPlotVector (kapa, xNew->n, xNew->data.F32, "x"); + KapaPlotVector (kapa, xNew->n, yNew->data.F32, "y"); + + pmVisualAskUser(&plotFixChips, &isVisual); + psFree(xNew); + psFree(yNew); + psFree(view); + return true; +} + + +bool pmAstromVisualPlotAstromGuessCheck (psVector *cornerPo, // P coordinates of chip corners before fitting + psVector *cornerQo, // Q coordinates of chip corners before fitting + psVector *cornerPn, // P coordinates of chip corners after fitting + psVector *cornerQn, // Q coordinates of chip corners after fitting + psVector *cornerPd, // P coordinate residuals of fit from old to new coordinates + psVector *cornerQd // Q coordinate residuals of fit from old to new coordinates + ) +{ + + //make sure we want to plot this + if (!isVisual || !plotAstromGuessCheck) return true; + + //set up graph window + if ( !pmVisualInitWindow (&kapa, "psastro:plots")) { + isVisual = false; + return false; + } + + Graphdata graphdata; + KapaSection section; + KapaInitGraph(&graphdata); + KapaClearPlots (kapa); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 2; + graphdata.size = 0.5; + graphdata.style = 2; + + section.dx = 0.4; + section.dy = 0.4; + + //Old Corners + section.x = 0.30; + section.y = 0.50; + section.name = NULL; + psStringAppend (§ion.name, "a0"); + KapaSetSection (kapa, §ion); + psFree(section.name); + + pmVisualScaleGraphdata (&graphdata, cornerPo, cornerPo, true); + KapaSetLimits (kapa, &graphdata); + KapaBox (kapa, &graphdata); + KapaSendLabel (kapa, "P (Pixels)", KAPA_LABEL_XM); + KapaSendLabel (kapa, "Q (Pixels)", KAPA_LABEL_YM); + KapaSendLabel (kapa, + "Fiducial Points in the Tangent Plane. Black: Initial Astrometry. Red: Final Astrometry", + KAPA_LABEL_XP); + KapaPrepPlot (kapa, cornerPo->n, &graphdata); + KapaPlotVector (kapa, cornerPo->n, cornerPo->data.F32, "x"); + KapaPlotVector (kapa, cornerQo->n, cornerQo->data.F32, "y"); + + // New Corners + graphdata.color = KapaColorByName("red"); + graphdata.ptype = 7; + graphdata.size = 1.5; + KapaPrepPlot (kapa, cornerPn->n, &graphdata); + KapaPlotVector (kapa, cornerPn->n, cornerPn->data.F32, "x"); + KapaPlotVector (kapa, cornerQn->n, cornerQn->data.F32, "y"); + + // Residuals + psVector *xResid = psVectorAlloc(cornerPn->n, PS_DATA_F32); + psVector *yResid = psVectorAlloc(cornerQn->n, PS_DATA_F32); + for(int i=0; i < cornerPn->n; i++) { + xResid->data.F32[i] = (cornerPd->data.F32[i]); + yResid->data.F32[i] = (cornerQd->data.F32[i]); + } + + graphdata.color = KapaColorByName("black"); + graphdata.size=0.5; + section.x = 0.3; + section.y = 0.0; + section.name = NULL; + psStringAppend (§ion.name, "a1"); + KapaSetSection (kapa, §ion); + psFree(section.name); + + pmVisualScaleGraphdata (&graphdata, xResid, yResid, true); + KapaSetLimits (kapa, &graphdata); + KapaBox (kapa, &graphdata); + KapaSendLabel (kapa, "dP", KAPA_LABEL_XM); + KapaSendLabel (kapa, "dQ", KAPA_LABEL_YM); + KapaSendLabel (kapa, + "Residual of the Fit from the Initial Astrometry to the Final Astrometry", + KAPA_LABEL_XP); + KapaPrepPlot (kapa, cornerPd->n, &graphdata); + KapaPlotVector (kapa, cornerPd->n, xResid->data.F32, "x"); + KapaPlotVector (kapa, cornerQd->n, yResid->data.F32, "y"); + + psFree(xResid); + psFree(yResid); + pmVisualAskUser (&plotAstromGuessCheck, &isVisual); + return true; +} + + +bool pmAstromVisualPlotCommonScale (pmFPA *fpa, // the fpa + psVector *oldScale // the old pixel scale of each chip in the fpa + ) +{ + //make sure we want to plot this + if (!isVisual || !plotCommonScale) return true; + + if (!pmVisualInitWindow(&kapa, "psastro:plots")){ + isVisual = false; + return false; + } + + KapaSection section = {"s1", .05, .05, .9, .9}; + Graphdata graphdata; + psPlane ptCH, ptFP; + ptCH.x = 0; + ptCH.y = 0; + psVector *xVec = psVectorAlloc (oldScale->n, PS_TYPE_F32); + psVector *yVec = psVectorAlloc (oldScale->n, PS_TYPE_F32); + + int nobj = 0; + + //project each chip corner to the Focal Plane + for(int i = 0; i < fpa->chips->n; i++) { + pmChip *chip = fpa->chips->data[i]; + if (!chip->process || !chip->file_exists) { continue; } + if (!chip->toFPA) { continue; } + + psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); + xVec->data.F32[nobj] = ptFP.x; + yVec->data.F32[nobj] = ptFP.y; + nobj++; + if (nobj == oldScale->n) break; + } + + //set up plot window + KapaInitGraph (&graphdata); + KapaClearPlots (kapa); + KapaSetSection (kapa, §ion); + KapaSetFont (kapa, "helvetica", 14); + pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, oldScale, false); + KapaSendLabel (kapa, "L (FP)", KAPA_LABEL_XM); + KapaSendLabel (kapa, "M (FP)", KAPA_LABEL_YM); + KapaSendLabel (kapa, "Old Pixel Scale of FPA Chips (Not to Scale)", KAPA_LABEL_XP); + + pmVisualAskUser (&plotCommonScale, &isVisual); + + psFree(xVec); + psFree(yVec); + + return true; +} + + +bool pmAstromVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars, + psArray *match, psMetadata *recipe) +{ + + //make sure we want to plot this + if (!isVisual || !plotMosaicOneChip) return true; + + if(!pmVisualInitWindow(&kapa, "psastro:plots") || !pmVisualInitWindow(&kapa2, "psastro:plots")) { + isVisual = false; + return false; + } + + //plot the residuals + if (!residPlot(rawstars, refstars, match, recipe, "Single Chip Fit Residuals - Mosaic Mode")) { + isVisual = false; + return false; + } + + //ask for user input and finish + pmVisualAskUser(&plotMosaicOneChip, &isVisual); + + return true; +} + + +bool pmAstromVisualPlotMosaicMatches( psArray *rawstars, psArray *refstars, + psArray *match, int iteration, + psMetadata *recipe) +{ + //make sure we want to plot this + if (!isVisual || !plotMosaicMatches) return true; + + char title[60]; + sprintf(title, "Matches found during psastroMosaicSetMatch iteration %d", iteration); + + + if(!pmVisualInitWindow(&kapa, "psastro:plots") || !pmVisualInitWindow(&kapa2, "psastro:plots")) { + isVisual = false; + return false; + } + + if (!residPlot(rawstars, refstars, match, recipe, title)){ + isVisual = false; + return false; + } + + //ask for user input + pmVisualAskUser (&plotMosaicMatches, &isVisual); + return true; +} + + +bool pmAstromVisualPlotGridMatch (const psArray *raw, + const psArray *ref, + psImage *gridNP, + double offsetX, + double offsetY, + double maxOffpix, + double Scale, + double Offset) { //make sure we want to plot this if (!isVisual || !plotGridMatch) return true; - if (!pmAstromVisualInitWindow(&kapa, "pmAstrom:plots")) - return false; + if (!pmVisualInitWindow(&kapa, "psastro:plots")){ + isVisual = false; + return false; + } KapaSection section = {"s1", 0.05, 0.05, .75, .75}; @@ -123,6 +819,7 @@ Graphdata graphdata; int nplot = raw->n * ref->n; // number of points to plot - float dXplot[nplot]; // x data points - float dYplot[nplot]; // y data points + psVector *dXplot = psVectorAlloc (nplot, PS_TYPE_F32); // x data points + psVector *dYplot = psVectorAlloc (nplot, PS_TYPE_F32); // y data points + pmAstromObj *ob1; pmAstromObj *ob2; // shortcuts to the data in raw and ref psU32 **NP = gridNP->data.U32; // shortcut to the gridNP data @@ -168,6 +865,6 @@ dX = ob1->FP->x - ob2->FP->x; dY = ob1->FP->y - ob2->FP->y; - dXplot[(i * ref->n) + j] = dX; - dYplot[(i * ref->n) + j] = dY; + dXplot->data.F32[(i * ref->n) + j] = dX; + dYplot->data.F32[(i * ref->n) + j] = dY; } } @@ -191,6 +888,6 @@ //Plot the offsets KapaPrepPlot(kapa, nplot, &graphdata); - KapaPlotVector (kapa, nplot, dXplot, "x"); - KapaPlotVector (kapa, nplot, dYplot, "y"); + KapaPlotVector (kapa, nplot, dXplot->data.F32, "x"); + KapaPlotVector (kapa, nplot, dYplot->data.F32, "y"); //Overplot bounding box, peak of distribution @@ -255,25 +952,21 @@ KapaPlotVector (kapa, 2, yslice, "y"); - pmAstromVisualAskUser(&plotGridMatch); + pmVisualAskUser(&plotGridMatch, &isVisual); + psFree(dXplot); + psFree(dYplot); return true; } // end of pmAstromVisualPlotGridMatch -/** Plot the refinements made within pmAstromGridTweak. - * After pmAstromGridMatch finds the best rotaion/scale/offset between raw and reference stars - * within a coarse grid of rotations/scales, pmAstromGridTweak computes a higher precision - * estimate of the offset. It computes the 2 point correlation function between raw and ref - * stars along horizontal and vertical cuts through the first-guess offset. It finds the peak - * of these two profiles and adjusts the offset accordingly. This procedure plots the profiles. - */ -bool pmAstromVisualPlotTweak (psVector *xHist, ///< Smoothed Horizontal cut through the histogram - psVector *yHist, ///< Smoothed Vertical cut throug the histogram - int xBin, ///< X Bin index of the histogram peak - int yBin ///< Y bin index of the histogram peak +bool pmAstromVisualPlotTweak (psVector *xHist, // Smoothed Horizontal cut through the histogram + psVector *yHist, // Smoothed Vertical cut throug the histogram + int xBin, // X Bin index of the histogram peak + int yBin // Y bin index of the histogram peak ) { //make sure we want to plot this if (!isVisual || !plotTweak) return true; - if (!pmAstromVisualInitWindow(&kapa, "pmAstrom:plots")) { + if (!pmVisualInitWindow(&kapa, "psastro:plots")) { + isVisual = false; return false; } @@ -299,5 +992,5 @@ // plot the X histogram - pmAstromVisualScaleGraphdata(&graphdata, xIndices, xHist, false); + pmVisualScaleGraphdata(&graphdata, xIndices, xHist, false); KapaSetSection(kapa, §ion1); KapaSetLimits (kapa, &graphdata); @@ -326,5 +1019,5 @@ //plot the Y histogram - pmAstromVisualScaleGraphdata(&graphdata, yIndices, yHist, false); + pmVisualScaleGraphdata(&graphdata, yIndices, yHist, false); KapaSetSection(kapa, §ion2); KapaSetLimits (kapa, &graphdata); @@ -358,5 +1051,5 @@ KAPA_LABEL_XP); - pmAstromVisualAskUser(&plotTweak); + pmVisualAskUser(&plotTweak, &isVisual); psFree(xIndices); @@ -366,84 +1059,300 @@ -/** pmAstromScaleGraphdata - * Scale the graphdata structure based on x and y coordinates. Use sigma clipping to - * prevent outliers from making te plot region too big. - */ -bool pmAstromVisualScaleGraphdata(Graphdata *graphdata, psVector *xVec, psVector *yVec, bool clip) { - - graphdata->xmin = +FLT_MAX; - graphdata->xmax = -FLT_MAX; - graphdata->ymin = +FLT_MAX; - graphdata->ymax = -FLT_MAX; - - //determine standard deviation of xVec and yVec - psStats *statsX = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - psStats *statsY = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - psVectorStats (statsX, xVec, NULL, NULL, 0); - psVectorStats (statsY, yVec, NULL, NULL, 0); - - float xhi = statsX->sampleMedian + 3 *statsX->sampleStdev; - float xlo = statsX->sampleMedian - 3 *statsX->sampleStdev; - float yhi = statsY->sampleMedian + 3 *statsY->sampleStdev; - float ylo = statsY->sampleMedian - 3 *statsY->sampleStdev; - - // don't sigma clip - if (!clip) { - xhi = +FLT_MAX; - xlo = -FLT_MAX; - yhi = +FLT_MAX; - ylo = -FLT_MAX; - } - - // abort if there is no good data - if (!isfinite(xhi) || !isfinite(xlo) || !isfinite(yhi) || !isfinite(ylo)) { - graphdata->xmin = -1; - graphdata->ymin = -1; - graphdata->xmax = 1; - graphdata->ymax = 1; - psFree(statsX); - psFree(statsY); - return false; - } - - for(int i = 0; i < xVec->n; i++) { - if (!isfinite(xVec->data.F32[i])) continue; - if (xVec->data.F32[i] > xhi || xVec->data.F32[i] < xlo) continue; - graphdata->xmin = PS_MIN (graphdata->xmin, xVec->data.F32[i]); - graphdata->xmax = PS_MAX (graphdata->xmax, xVec->data.F32[i]); - } - - for (int i = 0; i < yVec->n; i++) { - if (!isfinite(xVec->data.F32[i])) continue; - if (yVec->data.F32[i] > yhi || yVec->data.F32[i] < ylo) continue; - graphdata->ymin = PS_MIN (graphdata->ymin, yVec->data.F32[i]); - graphdata->ymax = PS_MAX (graphdata->ymax, yVec->data.F32[i]); - } - - // add a whitespace border - float range = graphdata->xmax - graphdata->xmin; - if (range == 0) range = 1; - graphdata->xmin -= .05 * range; - graphdata->xmax += .05 * range; - - range = graphdata->ymax - graphdata->ymin; - if (range == 0) range = 1; - graphdata->ymin -= .05 * range; - graphdata->ymax += .05 * range; - - psFree (statsX); - psFree (statsY); +bool residPlot (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe, + char *title) { + + + //initialize graph information + Graphdata graphdata; + KapaSection section; + + KapaInitGraph (&graphdata); + KapaClearPlots (kapa); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 7; + graphdata.size = 0.5; + graphdata.style = 2; + + section.dx = 0.4; + section.dy = 0.4; + + //initialize and populate the plotting vectors + bool status = false; + float iMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MIN"); + float iMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MAX"); + float rMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MIN"); + float rMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MAX"); + + psVector *xVec = psVectorAlloc (match->n, PS_TYPE_F32); + psVector *yVec = psVectorAlloc (match->n, PS_TYPE_F32); + psVector *zVec = psVectorAlloc (match->n, PS_TYPE_F32); + + // X vs dX + section.x = 0.0; + section.y = 0.5; + section.name = NULL; + psStringAppend (§ion.name, "a0"); + KapaSetSection (kapa, §ion); + psFree (section.name); + + int n = 0; + for (int i = 0; i < match->n; i++) { + pmAstromMatch *pair = match->data[i]; + pmAstromObj *raw = rawstars->data[pair->raw]; + pmAstromObj *ref = refstars->data[pair->ref]; + + if (!isfinite(raw->Mag)) continue; + if (raw->Mag < iMagMin) continue; + if (raw->Mag > iMagMax) continue; + if (ref->Mag < rMagMin) continue; + if (ref->Mag > rMagMax) continue; + + xVec->data.F32[n] = raw->chip->x; + yVec->data.F32[n] = raw->chip->x - ref->chip->x; + zVec->data.F32[n] = raw->Mag; + n++; + } + xVec->n = yVec->n = zVec->n = n; + + KapaSendLabel (kapa, "X", KAPA_LABEL_XM); + KapaSendLabel (kapa, "dX", KAPA_LABEL_YM); + pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false); + + // X vs dY + section.x = 0.5; + section.y = 0.5; + section.name = NULL; + psStringAppend (§ion.name, "a1"); + KapaSetSection (kapa, §ion); + psFree (section.name); + + n = 0; + for (int i = 0; i < match->n; i++) { + pmAstromMatch *pair = match->data[i]; + pmAstromObj *raw = rawstars->data[pair->raw]; + pmAstromObj *ref = refstars->data[pair->ref]; + + if (!isfinite(raw->Mag)) continue; + if (raw->Mag < iMagMin) continue; + if (raw->Mag > iMagMax) continue; + if (ref->Mag < rMagMin) continue; + if (ref->Mag > rMagMax) continue; + + xVec->data.F32[n] = raw->chip->x; + yVec->data.F32[n] = raw->chip->y - ref->chip->y; + zVec->data.F32[n] = raw->Mag; + n++; + } + xVec->n = yVec->n = zVec->n = n; + + KapaSendLabel (kapa, "X", KAPA_LABEL_XM); + KapaSendLabel (kapa, "dY", KAPA_LABEL_YM); + pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false); + + // Y vs dX + section.x = 0.0; + section.y = 0.0; + section.name = NULL; + psStringAppend (§ion.name, "a2"); + KapaSetSection (kapa, §ion); + psFree (section.name); + + n = 0; + for (int i = 0; i < match->n; i++) { + pmAstromMatch *pair = match->data[i]; + pmAstromObj *raw = rawstars->data[pair->raw]; + pmAstromObj *ref = refstars->data[pair->ref]; + + if (!isfinite(raw->Mag)) continue; + if (raw->Mag < iMagMin) continue; + if (raw->Mag > iMagMax) continue; + if (ref->Mag < rMagMin) continue; + if (ref->Mag > rMagMax) continue; + + xVec->data.F32[n] = raw->chip->y; + yVec->data.F32[n] = raw->chip->x - ref->chip->x; + zVec->data.F32[n] = raw->Mag; + n++; + } + xVec->n = yVec->n = zVec->n = n; + + KapaSendLabel (kapa, "Y", KAPA_LABEL_XM); + KapaSendLabel (kapa, "dX", KAPA_LABEL_YM); + pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false); + + // Y vs dY + section.x = 0.5; + section.y = 0.0; + section.name = NULL; + psStringAppend (§ion.name, "a3"); + KapaSetSection (kapa, §ion); + psFree (section.name); + + n = 0; + for (int i = 0; i < match->n; i++) { + pmAstromMatch *pair = match->data[i]; + pmAstromObj *raw = rawstars->data[pair->raw]; + pmAstromObj *ref = refstars->data[pair->ref]; + + if (!isfinite(raw->Mag)) continue; + if (raw->Mag < iMagMin) continue; + if (raw->Mag > iMagMax) continue; + if (ref->Mag < rMagMin) continue; + if (ref->Mag > rMagMax) continue; + + xVec->data.F32[n] = raw->chip->y; + yVec->data.F32[n] = raw->chip->y - ref->chip->y; + zVec->data.F32[n] = raw->Mag; + n++; + } + xVec->n = yVec->n = zVec->n = n; + + KapaSendLabel (kapa, "Y", KAPA_LABEL_XM); + KapaSendLabel (kapa, "dY", KAPA_LABEL_YM); + pmKapaPlotVectorTriple_AutoLimits_OpenGraph (kapa, &graphdata, xVec, yVec, zVec, false); + + section.x = 0.0; + section.y = 0.0; + section.dx = 0.95; + section.dy = 0.95; + section.name = NULL; + psStringAppend (§ion.name, "a5"); + KapaSetSection (kapa, §ion); + KapaSendLabel (kapa, title, KAPA_LABEL_XP); + psFree (section.name); + + //second window + + KapaInitGraph (&graphdata); + KapaClearPlots (kapa2); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 2; + graphdata.style = 2; + + psFree (xVec); + psFree (yVec); + psFree (zVec); + + xVec = psVectorAlloc (rawstars->n, PS_TYPE_F32); + yVec = psVectorAlloc (rawstars->n, PS_TYPE_F32); + zVec = psVectorAlloc (rawstars->n, PS_TYPE_F32); + + // X vs Y by mag (raw) + n = 0; + for (int i = 0; i < rawstars->n; i++) { + pmAstromObj *raw = rawstars->data[i]; + if (!isfinite(raw->Mag)) continue; + if (raw->Mag < iMagMin) continue; + if (raw->Mag > iMagMax) continue; + + xVec->data.F32[n] = raw->chip->x; + yVec->data.F32[n] = raw->chip->y; + zVec->data.F32[n] = raw->Mag; + n++; + } + xVec->n = yVec->n = zVec->n = n; + + KapaSendLabel (kapa2, "X", KAPA_LABEL_XM); + KapaSendLabel (kapa2, "Y", KAPA_LABEL_YM); + KapaSendLabel (kapa2, + "Chip Coordinates. Black = Raw Stars. Red = Ref Stars. Blue = Matched Stars" + , KAPA_LABEL_XP); + pmVisualTriplePlot (kapa2, &graphdata, xVec, yVec, zVec, false); + + // X vs Y by mag (ref) + psFree (xVec); + psFree (yVec); + psFree (zVec); + + xVec = psVectorAlloc (refstars->n, PS_TYPE_F32); + yVec = psVectorAlloc (refstars->n, PS_TYPE_F32); + zVec = psVectorAlloc (refstars->n, PS_TYPE_F32); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 7; + graphdata.style = 2; + + n = 0; + for (int i = 0; i < refstars->n; i++) { + pmAstromObj *ref = refstars->data[i]; + if (!isfinite(ref->Mag)) continue; + if (ref->Mag < rMagMin) continue; + if (ref->Mag > rMagMax) continue; + + xVec->data.F32[n] = ref->chip->x; + yVec->data.F32[n] = ref->chip->y; + zVec->data.F32[n] = ref->Mag; + n++; + } + xVec->n = yVec->n = zVec->n = n; + pmVisualTripleOverplot (kapa2, &graphdata, xVec, yVec, zVec, false); + + //rescale the graph to include all points + float xmin = graphdata.xmin; + float ymin = graphdata.ymin; + float xmax = graphdata.xmax; + float ymax = graphdata.ymax; + pmVisualScaleGraphdata(&graphdata, xVec, yVec, true); + graphdata.xmin = PS_MIN(xmin, graphdata.xmin); + graphdata.ymin = PS_MIN(ymin, graphdata.ymin); + graphdata.xmax = PS_MAX(xmax, graphdata.xmax); + graphdata.ymax = PS_MAX(ymax, graphdata.ymax); + KapaSetLimits (kapa2, &graphdata); + + //overplot matched stars in blue + psFree (xVec); + psFree (yVec); + psFree (zVec); + + xVec = psVectorAlloc (match->n, PS_TYPE_F32); + yVec = psVectorAlloc (match->n, PS_TYPE_F32); + zVec = psVectorAlloc (match->n, PS_TYPE_F32); + + graphdata.color = KapaColorByName ("blue"); + n = 0; + for (int i = 0; i < match->n; i++) { + pmAstromMatch *pair = match->data[i]; + pmAstromObj *raw = rawstars->data[pair->raw]; + pmAstromObj *ref = refstars->data[pair->ref]; + if (raw->Mag < iMagMin) continue; + if (raw->Mag > iMagMax) continue; + if (ref->Mag < rMagMin) continue; + if (ref->Mag > rMagMax) continue; + + xVec->data.F32[n] = raw->chip->x; + yVec->data.F32[n] = raw->chip->y; + zVec->data.F32[n] = iMagMin; + n++; + } + xVec->n = yVec->n = zVec->n = n; + pmVisualTripleOverplot (kapa2, &graphdata, xVec, yVec, zVec, false); + + psFree (xVec); + psFree (yVec); + psFree (zVec); return true; } + + # else bool pmAstromSetVisual(bool mode) { return true; } -bool pmAstromVisualInitWindow (int *kapid, char *name) { return true; } -bool pmAstromVisualAskUser (bool *plotflag) { return true; } bool pmAstromVisualClose() { return true; } bool pmAstromVisualPlotGridMatch (const psArray *raw, const psArray *ref, psImage *gridNP, double offsetX, double offsetY, double maxOffpix, double Scale, double Offset) { return true; } bool pmAstromVisualPlotTweak (psVector *xHist, psVector *yHist, int xBin, int yBin) {return true;} +bool pmAstromVisualPlotLuminosityFunction (psVector *lnMag, psVector *Mag, pmLumFunc *lumFunc, pmLumFunc *rawFunc) {return true;} +bool pmAstromVisualPlotRawStars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe) {return true;} +bool pmAstromVisualPlotRefStars (psArray *refstars, psMetadata *recipe) {return true;} +bool pmAstromVisualPlotRemoveClumps (psArray *input, psImage *count, int scale, float limit) {return true;} +bool pmAstromVisualPlotFixChips (pmFPAfile *input, psVector *xOld, psVector *yOld) {return true;} +bool pmAstromVisualPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe) {return true;} +bool pmAstromVisualPlotAstromGuessCheck (psVector *cornerPo, psVector *cornerQo, psVector *cornerPn, psVector *cornerQn, psVector *cornerPd, psVector *cornerQd) {return true;} +bool pmAstromVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe) {return true;} +bool pmAstromVisualPlotCommonScale (pmFPA *fpa, psVector *oldScale) {return true;} +bool pmAstromVisualPlotMosaicMatches (psArray *rawstars, psArray *refstars, psArray *match, int iteration, psMetadata *recipe) {return true;} # endif Index: trunk/psModules/src/astrom/pmAstrometryVisual.h =================================================================== --- trunk/psModules/src/astrom/pmAstrometryVisual.h (revision 21387) +++ trunk/psModules/src/astrom/pmAstrometryVisual.h (revision 21422) @@ -1,6 +1,169 @@ +/* + * @file pmAstrometryVisual.h + * @author Chris Beaumont, IfA + * @brief A set of functions to display visual diagnostics from psastro + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_ASTROM_VISUAL_H +#define PM_ASTROM_VISUAL_H + + +/** A fit to the logN / logS curve for a set of stars + * logN = offset + slope * logS + */ +typedef struct { + double mMin; ///< minimum magnitude bin with data + double mMax; ///< maximum magnitude bin with data + double offset; ///< fitted line offset + double slope; ///< fitted line slope + double mPeak; ///< mag of peak bin + int nPeak; ///< # of stars in peak bin + int sPeak; ///< sum of stars to peak bin +} pmLumFunc; + + +/** Enable or disable visual plotting for psastro routines + * @param mode true/false to enable/disable plotting + * @return true for success */ bool pmAstromSetVisual(bool mode); -bool pmAstromVisualInitWindow (int *kapid, char *name); -bool pmAstromVisualAskUser (bool *plotflag); + + +/** Close plotting windows at the end of a run + * @return true for success */ bool pmAstromVisualClose(); -bool pmAstromVisualPlotGridMatch (const psArray *raw, const psArray *ref, psImage *gridNP, double offsetX, double offsetY, double maxOffpix, double Scale, double Offset); -bool pmAstromVisualPlotTweak (psVector *xHist, psVector *yHist, int xBin, int yBin); + + +/** + * Plot the offset between every pair of reference and raw source locations. The peak of this + * distribution nominally gives the offset, scale difference, and rotation of the two catalogs. + * Overplots the location of this peak as determined by pmAstromGridAngle, as well as some profiles + * along horizontal and vertical cuts through this peak. + */ +bool pmAstromVisualPlotGridMatch (const psArray *raw, ///< raw stars + const psArray *ref, ///< reference stars + psImage *gridNP, ///< a 2D histogram of raw-ref star distances + double offsetX, ///< The X location (FP coordinates) of the peak of gridNP + double offsetY, ///< the Y location (FP coordinates) of the peak of gridNP + double maxOffpix, ///< The half-width of gridNP in FP coordinates + double Scale, ///< The pixel size of gridNP in histogram-bin-coordinates + double Offset ///< The (x,y) location (histogram-bin coordinates) of the FP point (0,0) in gridNP + ); + + +/** + * Plot the refinements made within pmAstromGridTweak. + * After pmAstromGridMatch finds the best rotaion/scale/offset between raw and reference stars + * within a coarse grid of rotations/scales, pmAstromGridTweak computes a higher precision + * estimate of the offset. It computes the 2 point correlation function between raw and ref + * stars along horizontal and vertical cuts through the first-guess offset. It finds the peak + * of these two profiles and adjusts the offset accordingly. This procedure plots the profiles. + */ +bool pmAstromVisualPlotTweak (psVector *xHist, ///< Smoothed Horizontal cut through the histogram + psVector *yHist, ///< Smoothed Vertical cut throug the histogram + int xBin, ///< X Bin index of the histogram peak + int yBin ///< Y bin index of the histogram peak + ); + + +/** + * Plot the two luminosity functions created within psastroRefStarSubset + * The luminosity functions are used to select a subset of reference stars, + * so we plot the cutoff that defines this subset + */ +bool pmAstromVisualPlotLuminosityFunction (psVector *lnMag, ///< Log(n) for each magnitude bin + psVector *Mag, ///< magnitude bins + pmLumFunc *lumFunc,///< Fit to the reference star luminosity function + pmLumFunc *rawFunc ///< Fit to the raw star luminoisty function + ); + + +/** + * Plot raw stars as determined from first pass astrometry fit + * Called within psastroAstromGeuss + */ +bool pmAstromVisualPlotRawStars (psArray *rawstars, ///< Stars detected in the fpa + pmFPA *fpa, ///< structure describing the focal plane array + pmChip *chip, ///< structure describing the chip + psMetadata *recipe ///< the recipe used in psastro + ); + + +/** + * plot the location of references stars over the entire fpa + * invoked during psastroChooseRefStars + */ +bool pmAstromVisualPlotRefStars (psArray *refstars, psMetadata *recipe); + + +/** + * Plot the stars in a region, and indicate which stars are part of 'clumps' + * These stars are flagged during astrometric fitting, since dense regions are + * harder to cross-match than sparse ones. Called during psastroRemoveClumps. + */ +bool pmAstromVisualPlotRemoveClumps (psArray *input, ///< Array containing the field stars + psImage *count, ///< A 2D histogram of the field star distribution + int scale, ///< The pixel size of the histogram + float limit ///< The minimum numuber of stars in a bin flagged as a clump + ); + +/** + * Plots the chip corners in the FP before and after chips with inconsistent solutions have been fixed. + * Invoked during psastroFixChips + */ +bool pmAstromVisualPlotFixChips (pmFPAfile *input, ///< focal plane array file + psVector *xOld, ///< old X location of chip cornerss + psVector *yOld ///< old Y location of chip corners + ); + + +/** + * Assess the goodness of fit for a signle chip by + * plotting the fit residuals + * invoked during psastroOneChipFit + */ +bool pmAstromVisualPlotOneChipFit (psArray *rawstars, ///< stars detected in the image + psArray *refstars, ///< reference stars over the same region + psArray *match, ///< contains which rawstars match to which refstars + psMetadata *recipe ///< data reduction recipe + ); + +/** + * Plots the fpa chip corners projected on to the tangential plane before and after + * the astrometry solution has been applied. In psastroAstromGuessCheck, the old corners + * are then fit to the new corners to get a sense at how far off the initial WCS info was + * in offset, rotation, and scale. This procedure also plots the residuals of the fit from + * old to new coordinates + */ +bool pmAstromVisualPlotAstromGuessCheck (psVector *cornerPo, ///< P coordinates of chip corners before fitting + psVector *cornerQo, ///< Q coordinates of chip corners before fitting + psVector *cornerPn, ///< P coordinates of chip corners after fitting + psVector *cornerQn, ///< Q coordinates of chip corners after fitting + psVector *cornerPd, ///< P coordinate residuals of fit from old to new coordinates + psVector *cornerQd ///< Q coordinate residuals of fit from old to new coordinates + ); + + +/** + * plot the residuals between raw stars and ref stars after + * fitting in psastroMosaicOneChip + */ +bool pmAstromVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe) ; + + +/** + * Plots the pixel scales of the fpa before they are + * equalized in psastroMosaicCommonScale + */ +bool pmAstromVisualPlotCommonScale (pmFPA *fpa, ///< the fpa + psVector *oldScale ///< the old pixel scale of each chip in the fpa + ); + + +/** pmAstromVisualPlotMosaicMatches + * Plot the matches between raw and reference stars during pmAstromVisualMosaicSetMatch + */ +bool pmAstromVisualPlotMosaicMatches (psArray *rawstars, psArray *refstars, psArray *match, int iteration, psMetadata *recipe); + + +#endif Index: trunk/psModules/src/extras/.cvsignore =================================================================== --- trunk/psModules/src/extras/.cvsignore (revision 21387) +++ trunk/psModules/src/extras/.cvsignore (revision 21422) @@ -3,7 +3,4 @@ Makefile Makefile.in -libpsmodulesextras.la -libpsmodulesextras_la-pmKapaPlots.lo -libpsmodulesextras_la-psIOBuffer.lo -libpsmodulesextras_la-psPipe.lo -libpsmodulesextras_la-psVectorBracket.lo +*.la +*.lo Index: trunk/psModules/src/extras/Makefile.am =================================================================== --- trunk/psModules/src/extras/Makefile.am (revision 21387) +++ trunk/psModules/src/extras/Makefile.am (revision 21422) @@ -7,5 +7,6 @@ psPipe.c \ psIOBuffer.c \ - pmKapaPlots.c + pmKapaPlots.c \ + pmVisual.c pkginclude_HEADERS = \ @@ -13,5 +14,6 @@ psPipe.h \ psIOBuffer.h \ - pmKapaPlots.h + pmKapaPlots.h \ + pmVisual.h CLEANFILES = *~ Index: trunk/psModules/src/extras/pmVisual.c =================================================================== --- trunk/psModules/src/extras/pmVisual.c (revision 21422) +++ trunk/psModules/src/extras/pmVisual.c (revision 21422) @@ -0,0 +1,327 @@ +/** The following are a collection of core procedures to aid the creation of visaual diagnostics + * @author Chris Beaumont, IfA + * @date January 23, 2008 + **/ + +/* Include Files */ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmAstrometryObjects.h" +#include "pmFPAExtent.h" + +# if (HAVE_KAPA) +# include +# include "pmVisual.h" +# include "pmKapaPlots.h" +# define KAPAX 700 +# define KAPAY 700 + + +bool pmVisualInitWindow (int *kapid, char *name) { + if (*kapid == -1) { + *kapid = KapaOpenNamedSocket("kapa", name); + if (*kapid == -1) { + fprintf (stderr, "Failure to open kapa.\n"); + return false; + } + KapaResize (*kapid, KAPAX, KAPAY); + } + return true; +} + + +bool pmVisualInitGraph (int kapa, KapaSection *section, Graphdata *graphdata) +{ + KapaSetSection (kapa, section); + KapaSetFont (kapa, "helvetica", 14); + KapaSetLimits (kapa, graphdata); + KapaBox (kapa, graphdata); + return true; +} + + +bool pmVisualAskUser(bool *plotFlag, bool *packageFlag) +{ + char key[10]; + fprintf (stdout, "[c]ontinue? [s]kip the rest of these plots? [a]bort all visual plots?"); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + if (key[0] == 's') { + *plotFlag = false; + } + if (key[0] == 'a') { + *packageFlag = false; + } + return true; +} + + +bool pmVisualImStats(psImage *image, double *mean, double *stdev, double *min, double *max) { + + *min = +FLT_MAX; + *max = -FLT_MAX; + double ex = 0; // + double ex2 = 0; // + int numPix = 0; // number of finite pixels + bool isDouble = image->type.type == PS_TYPE_F64; + + if(!isDouble && (image->type.type != PS_TYPE_F32)) { + fprintf(stderr, "Image must be PS_TYPE_F32 or PS_TYPE_F64\n"); + return false; + } + + for(int i = 0; i < image->numRows; i++) { + for(int j = 0; j < image->numCols; j++) { + double entry; + if(isDouble) { + if (!isfinite(image->data.F64[i][j])) continue; + entry = image->data.F64[i][j]; + } else { + if (!isfinite(image->data.F32[i][j])) continue; + entry = image->data.F32[i][j]; + } + numPix++; + ex += entry; + ex2 += (entry * entry); + *min = PS_MIN(*min, entry); + *max = PS_MAX(*max, entry); + } + } + + if (numPix == 0) return false; + ex /= numPix; + ex2 /= numPix; + *mean = ex; + *stdev = sqrt(ex2 - ex * ex); + + return true; +} + + +bool pmVisualTriplePlot (int kapid, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing) +{ + pmVisualScaleGraphdata (graphdata, xVec, yVec, true); + //printf("%f %f %f %f \n",graphdata->xmin, graphdata->xmax, graphdata->ymin, graphdata->ymax); + // set the scale vector + psVector *zScale = psVectorAlloc (zVec->n, PS_DATA_F32); + pmVisualCreateScaleVec (zVec, zScale, increasing); + + KapaSetFont (kapid, "helvetica", 14); + KapaSetLimits(kapid, graphdata); + KapaBox (kapid, graphdata); + + // the point size will be scaled from the z vector + graphdata->ptype = 7; + graphdata->style = 2; + graphdata->size = -1; + KapaPrepPlot (kapid, xVec->n, graphdata); + KapaPlotVector (kapid, xVec->n, xVec->data.F32, "x"); + KapaPlotVector (kapid, yVec->n, yVec->data.F32, "y"); + KapaPlotVector (kapid, zVec->n, zScale->data.F32, "z"); + psFree (zScale); + return true; +} + + +bool pmVisualTripleOverplot (int kapid, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing) { + + // set the scale vector + psVector *zScale = psVectorAlloc (zVec->n, PS_DATA_F32); + pmVisualCreateScaleVec (zVec, zScale, increasing); + + KapaSetFont (kapid, "helvetica", 14); + + // the point size will be scaled from the z vector + graphdata->size = -1; + KapaPrepPlot (kapid, xVec->n, graphdata); + KapaPlotVector (kapid, xVec->n, xVec->data.F32, "x"); + KapaPlotVector (kapid, yVec->n, yVec->data.F32, "y"); + KapaPlotVector (kapid, zVec->n, zScale->data.F32, "z"); + psFree (zScale); + return true; +} + + +bool pmVisualCreateScaleVec (psVector *zVec, psVector *zScale, bool increasing) { + // set limits based on data values + float zmin = +FLT_MAX; + float zmax = -FLT_MAX; + + for (int i = 0; i < zVec->n; i++) { + zmin = PS_MIN (zmin, zVec->data.F32[i]); + zmax = PS_MAX (zmax, zVec->data.F32[i]); + } + + float range = zmax - zmin; + if (range == 0.0) { + psVectorInit (zScale, 1.0); + } else { + for (int i = 0; i < zVec->n; i++) { + if (increasing) { + zScale->data.F32[i] = PS_MIN (1.5, PS_MAX(0.05, 1.5*(zVec->data.F32[i] - zmin)/range)); + } else { + zScale->data.F32[i] = PS_MIN (1.5, PS_MAX(0.05, 1.5*(zmax - zVec->data.F32[i])/range)); + } + } + } + return true; +} + + +bool pmVisualScaleImage(int kapaFD, psImage *inImage, const char *name, int channel, bool clip) { + + KiiImage image; + KapaImageData data; + Coords coords; + + //make sure we have a compatible image type + if(inImage->type.type != PS_TYPE_F32) { + fprintf(stderr, "Cannot display this image (imcompatible data type)\n"); + return false; + } + + strcpy (coords.ctype, "RA---TAN"); + + double min, max, stdev, mean; + if(!pmVisualImStats(inImage, &mean, &stdev, &min, &max)) 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 = clip ? mean - 3 * stdev : min; + data.range = clip ? 6 * stdev : max - min; + data.logflux = 0; + + KiiSetChannel (kapaFD, channel); + KiiNewPicture2D (kapaFD, &image, &data, &coords); + + return true; +} + + +psImage* pmVisualImageToFloat(psImage *image) { + psImage *result = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32); + + if (image->type.type == PS_TYPE_F32) { + psImageOverlaySection(result, image, 0, 0, "="); + return image; + } else if (image->type.type == PS_TYPE_F64) { + for (int i = 0; i < image->numRows; i++) { + for (int j = 0; j < image->numCols; j++) { + result->data.F32[i][j] = (float) image->data.F64[i][j]; + } + } + } else { + fprintf(stderr, "Unsupported psImage data type (F32 or F64)"); + return NULL; + } + return result; +} + + +bool pmVisualScaleGraphdata(Graphdata *graphdata, psVector *xVec, psVector *yVec, bool clip) { + + graphdata->xmin = +FLT_MAX; + graphdata->xmax = -FLT_MAX; + graphdata->ymin = +FLT_MAX; + graphdata->ymax = -FLT_MAX; + + //determine standard deviation of xVec and yVec + psStats *statsX = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + psStats *statsY = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + psVectorStats (statsX, xVec, NULL, NULL, 0); + psVectorStats (statsY, yVec, NULL, NULL, 0); + + float xhi = statsX->sampleMedian + 3 *statsX->sampleStdev; + float xlo = statsX->sampleMedian - 3 *statsX->sampleStdev; + float yhi = statsY->sampleMedian + 3 *statsY->sampleStdev; + float ylo = statsY->sampleMedian - 3 *statsY->sampleStdev; + + // don't sigma clip + if (!clip) { + xhi = +FLT_MAX; + xlo = -FLT_MAX; + yhi = +FLT_MAX; + ylo = -FLT_MAX; + } + + for(int i = 0; i < xVec->n; i++) { + if (!isfinite(xVec->data.F32[i])) continue; + if (xVec->data.F32[i] > xhi || xVec->data.F32[i] < xlo) continue; + graphdata->xmin = PS_MIN (graphdata->xmin, xVec->data.F32[i]); + graphdata->xmax = PS_MAX (graphdata->xmax, xVec->data.F32[i]); + } + + for (int i = 0; i < yVec->n; i++) { + if (!isfinite(xVec->data.F32[i])) continue; + if (yVec->data.F32[i] > yhi || yVec->data.F32[i] < ylo) continue; + graphdata->ymin = PS_MIN (graphdata->ymin, yVec->data.F32[i]); + graphdata->ymax = PS_MAX (graphdata->ymax, yVec->data.F32[i]); + } + + + // abort if there is no good data + if (!isfinite(xhi) || !isfinite(xlo) || !isfinite(yhi) || !isfinite(ylo)) { + graphdata->xmin = -1; + graphdata->ymin = -1; + graphdata->xmax = 1; + graphdata->ymax = 1; + psFree(statsX); + psFree(statsY); + return false; + } + + // add a whitespace border + float range = graphdata->xmax - graphdata->xmin; + if (range == 0) range = 1; + graphdata->xmin -= .05 * range; + graphdata->xmax += .05 * range; + + range = graphdata->ymax - graphdata->ymin; + if (range == 0) range = 1; + graphdata->ymin -= .05 * range; + graphdata->ymax += .05 * range; + + psFree (statsX); + psFree (statsY); + return true; +} + +#else + +bool pmVisualInitWindow(int *kapid, char *name){return true;} +bool pmVisualInitGraph (int kapa, void *section, void *graphdata){return true;} + +bool pmVisualAskUser(bool *plotFlag, bool *packageFlag){return true;} +bool pmVisualImStats(psImage *image, double *mean, + double *stdev, double *min, double *max){return true;} +bool pmVisualTriplePlot (int kapid, void *graphdata, psVector *xVec, + psVector *yVec, psVector *zVec, bool increasing){return true;} +bool pmVisualTripleOverplot (int kapid, void *graphdata, psVector *xVec, + psVector *yVec, psVector *zVec, bool increasing){return true;} +bool pmVisualCreateScaleVec (psVector *zVec, psVector *zScale, bool increasing){return true;} +bool pmVisualResidPlot (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe, char *title, int *kapa, int *kapa2){return true;} +bool pmVisualScaleImage(int kapaFD, psImage *inImage, + const char *name, int channel, bool clip){return true;} +bool pmVisualScaleGraphdata(void *graphdata, psVector *xVec, + psVector *yVec, bool clip){return true;} + +psImage* pmVisualImageToFloat(psImage *image){return NULL;} + + +#endif Index: trunk/psModules/src/extras/pmVisual.h =================================================================== --- trunk/psModules/src/extras/pmVisual.h (revision 21422) +++ trunk/psModules/src/extras/pmVisual.h (revision 21422) @@ -0,0 +1,137 @@ +/* @file pmVisual.h + * @brief functions to create visual diagnostics with the help of 'kapa' + * @author Chris Beaumont, IfA + * + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_VISUAL_H +#define PM_VISUAL_H + +#if (HAVE_KAPA) + + +/** Open, name, and resize a window for plotting. + * @param kapid an identifier for this window. Initialize to -1. Its value will be updated. + * @param name What to name the window. Seems not to like spaces. + * @return true for successful completion +*/ +bool pmVisualInitWindow (int *kapid, char *name); + + +/** Initialize a graph, set the appropriate section and font. + * @param kapa the index of the kapa window to plot to + * @param section the section to use + * @graphdata to use + * @return true for successful completion + */ +bool pmVisualInitGraph (int kapa, KapaSection *section, Graphdata *graphdata); + + +/** Ask the user how to proceed. + * At the user's request, this will disable diagnostic plotting. + * @param plotFlag, set to false if this plot should be disabled in the future + * @param packageFlag, set to false if all plots from this package (e.g. psastro, pmSubtraction) + * should be disabled. + */ +bool pmVisualAskUser(bool *plotFlag, bool *packageFlag); + + +/** Scale and display an image. + * @param kapaFD the index of the Kapa window to draw to + * @param inImage the image to display + * @param name the image label + * @param which channel to draw to (zero works) + * @param clip set to true to sigma clip the image when scaling + * @return true for successful completion */ +bool pmVisualScaleImage(int kapaFD, psImage *inImage, + const char *name, int channel, bool clip); + + +/** Calculate statistics on an image. + * This can handle non-finite input pixels + * @param image to calculate statistics for + * @param mean stores the calculated mean + * @param stdev stores the calculated standard deviation + * @param min stores the calculated minimum + * @param max stores the calculated maximum + * @return true for successful completion */ +bool pmVisualImStats(psImage *image, double *mean, + double *stdev, double *min, double *max); + + +/** Create a PS_TYPE_F32 image out of a PS_TYPE_F32 or PS_TYPE_F64 image. + * This is needed when displaying PS_TYPE_F64 images, which will not work wtih Kii. + * @param image to convert to type PS_TYPE_F32 + * @return a copy of the image of type PS_TYPE_F32, or NULL upon failure */ +psImage* pmVisualImageToFloat(psImage *image); + + +/** Create a scaled vector of plot point sizes from an unscaled, raw vector. + * Used for input into pmVisualTriplePlot + * @param zVec the raw data + * @param zScale the scaled vector of plot point sizes + * @return true for successful completion + */ +bool pmVisualCreateScaleVec (psVector *zVec, psVector *zScale, bool increasing); + + +/** Use x and y data to determine appropriate values for a Graphdata structure. + * This procedure sets the max and min keywords of a Graphdata structure to encompass + * the data coordinates given in xVec and yVec. Optionally, it will try to ignore outliers. + * @param graphdata structure to set up + * @param xVec X coordinates of data points + * @param yVec Y coordinates of data points + * @clip set to true to attempt to clip out outliers + * @return true for successful completion + */ +bool pmVisualScaleGraphdata(Graphdata *graphdata, psVector *xVec, + psVector *yVec, bool clip); + + +/** Create a scatter plot form a set of (x,y) positions, whose plot points + * are scaled by the size of a z vector + * @param kapid the index of the kapa window to plot to + * @param graphdata describing the plotting window + * @param xVec x positions + * @param yVec y positions + * @param zVec plot point sizes + * @param increasing whether the points in x,y,zvec are sorted + */ +bool pmVisualTriplePlot (int kapid, Graphdata *graphdata, psVector *xVec, + psVector *yVec, psVector *zVec, bool increasing); + +/** The same as pmVisualTriplePlot, but draws new points without erasing the old plot window + * @param kapid the index of the kapa window to plot to + * @param xVec x positions + * @param yVec y positions + * @param zVec plot point sizes + * @param increasing whether the x,y,zvectors are listed in increasing order + */ +bool pmVisualTripleOverplot (int kapid, Graphdata *graphdata, psVector *xVec, + psVector *yVec, psVector *zVec, bool increasing); + +#else + +// kapa-specific data types are changed to void +bool pmVisualInitWindow (int *kapid, char *name); +bool pmVisualInitGraph (int kapa, void *section, void *graphdata); +bool pmVisualAskUser(bool *plotFlag, bool *packageFlag); +bool pmVisualScaleImage(int kapaFD, psImage *inImage, + const char *name, int channel, bool clip); +bool pmVisualImStats(psImage *image, double *mean, + double *stdev, double *min, double *max); +psImage* pmVisualImageToFloat(psImage *image); +bool pmVisualCreateScaleVec (psVector *zVec, psVector *zScale, bool increasing); +bool pmVisualResidPlot (psArray *rawstars, psArray *refstars, psArray *match, + psMetadata *recipe, char *title, int *kapa, int *kapa2); +bool pmVisualScaleGraphdata(void *graphdata, psVector *xVec, + psVector *yVec, bool clip); +bool pmVisualTriplePlot (int kapid, void *graphdata, psVector *xVec, + psVector *yVec, psVector *zVec, bool increasing); +bool pmVisualTripleOverplot (int kapid, void *graphdata, psVector *xVec, + psVector *yVec, psVector *zVec, bool increasing); + +#endif //HAVE_KAPA + +#endif //ndef PM_VISUAL_H Index: trunk/psModules/src/imcombine/Makefile.am =================================================================== --- trunk/psModules/src/imcombine/Makefile.am (revision 21387) +++ trunk/psModules/src/imcombine/Makefile.am (revision 21422) @@ -18,5 +18,6 @@ pmSubtractionStamps.c \ pmSubtractionThreads.c \ - pmPSFEnvelope.c + pmPSFEnvelope.c \ + pmSubtractionVisual.c pkginclude_HEADERS = \ @@ -35,5 +36,6 @@ pmSubtractionStamps.h \ pmSubtractionThreads.h \ - pmPSFEnvelope.h + pmPSFEnvelope.h \ + pmSubtractionVisual.h CLEANFILES = *~ Index: trunk/psModules/src/imcombine/pmSubtractionAnalysis.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 21387) +++ trunk/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 21422) @@ -12,4 +12,5 @@ #include "pmSubtractionAnalysis.h" +#include "pmSubtractionVisual.h" #define KERNEL_MOSAIC 2 // Half-number of kernel instances in the mosaic image @@ -101,4 +102,5 @@ } + pmSubtractionVisualPlotConvKernels(convKernels); psMetadataAddImage(analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL.IMAGE", PS_META_DUPLICATE_OK, "Realisations of kernel", convKernels); Index: trunk/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 21387) +++ trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 21422) @@ -13,5 +13,5 @@ #include "pmSubtractionEquation.h" - +#include "pmSubtractionVisual.h" //#define TESTING @@ -695,4 +695,6 @@ } + pmSubtractionVisualPlotLeastSquares(stamps); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate equation: %f sec", psTimerClear("pmSubtractionCalculateEquation")); @@ -1001,4 +1003,5 @@ } + pmSubtractionVisualPlotLeastSquares((pmSubtractionStampList *) stamps); //casting away const return true; } Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 21387) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 21422) @@ -19,5 +19,5 @@ #include "pmSubtractionThreads.h" #include "pmSubtractionMatch.h" - +#include "pmSubtractionVisual.h" #define BG_STAT PS_STAT_ROBUST_MEDIAN // Statistic to use for background @@ -86,5 +86,5 @@ memCheck(" extract stamps"); - + pmSubtractionVisualPlotStamps(*stamps, (pmReadout *) ro1); return true; } Index: trunk/psModules/src/imcombine/pmSubtractionVisual.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionVisual.c (revision 21422) +++ trunk/psModules/src/imcombine/pmSubtractionVisual.c (revision 21422) @@ -0,0 +1,275 @@ +/** Diagnostic plots for pmSubtraction + * @author Chris Beaumont, IfA + */ + +/* Include Files */ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include + +#include "pmKapaPlots.h" +#include "pmSubtractionStamps.h" + +#include "pmVisual.h" + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmAstrometryObjects.h" + +# if (HAVE_KAPA) +# include + +//variables to determine when things are plotted +static bool isVisual = false; +static bool plotConvKernels = true; +static bool plotStamps = true; +static bool plotLeastSquares = true; +static bool plotImage = true; +// variables to store plotting window indices +static int kapa = -1; +static int kapa2 = -1; + +/** function prototypes*/ +static bool plotStampLocations(pmSubtractionStampList *stamps, pmReadout *ro); + +// Initialization Routines + +/** start or stop plotting */ +bool pmSubtractionSetVisual (bool mode) { + isVisual = mode; + return true; +} + + +/** destroy windows at the end of a run*/ +bool pmSubtractionVisualClose() +{ + if(kapa != -1) + KiiClose(kapa); + if(kapa2 != -1) + KiiClose(kapa2); + return true; +} + +// Plotting Routines + +/** Display images of the convolution kernels + * @param convKernels the kernels to plot + * @return true for success */ +bool pmSubtractionVisualPlotConvKernels(psImage *convKernels) { + if (!isVisual || !plotConvKernels) return true; + if (!pmVisualInitWindow(&kapa, "ppSub:Images")) { + isVisual = false; + return false; + } + pmVisualScaleImage(kapa, convKernels, "Convolution_Kernels", 0, true); + pmVisualAskUser(&plotConvKernels, &isVisual); + return true; +} + + +/** Display the postage stamps used to determine the convolution kernels + @param stamps to display + @return true for success */ +bool pmSubtractionVisualPlotStamps(pmSubtractionStampList *stamps, pmReadout *ro) { + if (!isVisual || !plotStamps) return true; + if (!pmVisualInitWindow (&kapa, "ppSub:Images")) { + isVisual = false; + return false; + } + if (!pmVisualInitWindow (&kapa2, "ppSub:StampMasterImage")) { + isVisual = false; + return false; + } + + //Plot Location of stamps + plotStampLocations(stamps, ro); + + //Find the stamp size + int imageMax = -1; + int numStamps = 0; + psElemType type = PS_TYPE_F32; // will be overwritten + for(int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; + if (stamp == NULL) continue; + + psKernel *k1 = stamp->image1; + if (k1 == NULL) continue; + + psImage *i1 = k1->image; + if (i1 == NULL) continue; + + imageMax = PS_MAX(imageMax, i1->numCols); + imageMax = PS_MAX(imageMax, i1->numRows); + numStamps++; + type = i1->type.type; + } + if (imageMax == -1) return false; + + int border = 30; + int tileRowCount = (int) ceil(sqrt(numStamps)); + int canvasX = tileRowCount * (imageMax * 2 + border); + int canvasY = tileRowCount * (imageMax + border); + psImage *canvas = psImageAlloc (canvasX, canvasY, type); + psImageInit (canvas, NAN); + + //overlay the images + int stampNum = 0; + int stampListNum = 0; + while (stampNum < numStamps) { + int x0 = (2 * imageMax + border) * (stampNum % tileRowCount); + int y0 = (imageMax + border) * (stampNum / tileRowCount); + + pmSubtractionStamp *stamp = stamps->stamps->data[stampListNum++]; + if (stamp == NULL) continue; + + psKernel *k1 = stamp->image1; + psKernel *k2 = stamp->image2; + + if (k1 == NULL) continue; + psImage *i1 = k1->image; + psImage *i2 = k2->image; + + if (i1 == NULL) continue; + psImageOverlaySection(canvas, i1, x0 + 0 * imageMax, y0, "="); + psImageOverlaySection(canvas, i2, x0 + 1 * imageMax + border / 4, y0, "="); + stampNum++; + } + pmVisualScaleImage(kapa, canvas, "Subtraction_Stamps", 0, true); + + pmVisualAskUser(&plotStamps, &isVisual); + return true; +} + +/** Plot the least-squares matrix of each stamp */ +bool pmSubtractionVisualPlotLeastSquares (pmSubtractionStampList *stamps, bool dual) { + + if (!isVisual || !plotLeastSquares) return true; + if (!pmVisualInitWindow (&kapa, "PPSub:Images")) { + isVisual = false; + return false; + } + + //Find the stamp size + int imageMax = -1; + int numStamps = 0; + psElemType type = PS_TYPE_F64; + + for(int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; + if (stamp == NULL) continue; + + psImage *im = stamp->matrix1; + if (im == NULL) im = stamp->matrix2; + if (im == NULL) im = stamp->matrixX; + if (im == NULL) continue; + + imageMax = PS_MAX(imageMax, im->numCols); + imageMax = PS_MAX(imageMax, im->numRows); + numStamps++; + type = im->type.type; + } + if (imageMax == -1) return false; + + int border = 15; + imageMax += border; + int tileRowCount = (int) ceil(sqrt(numStamps)); + int canvasX = tileRowCount * (imageMax); + int canvasY = tileRowCount * (imageMax); + psImage *canvas = psImageAlloc (canvasX, canvasY, type); + psImageInit (canvas, NAN); + + //overlay the images + int stampNum = 0; + int stampListNum = 0; + while (stampNum < numStamps) { + int x0 = (imageMax) * (stampNum % tileRowCount); + int y0 = (imageMax) * (stampNum / tileRowCount); + + pmSubtractionStamp *stamp = stamps->stamps->data[stampListNum++]; + if (stamp == NULL) continue; + + psImage *im = stamp->matrix1; + if (im == NULL) im = stamp->matrix2; + if (im == NULL) im = stamp->matrixX; + if (im == NULL) continue; + + psImageOverlaySection(canvas, im, x0, y0, "="); + stampNum++; + } + + psImage *canvas32 = pmVisualImageToFloat(canvas); + pmVisualScaleImage(kapa, canvas32, "Least_Squares", 0, true); + + pmVisualAskUser(&plotLeastSquares, &isVisual); + psFree(canvas); + psFree(canvas32); + return true; +} + +bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub) { + if (!isVisual || !plotImage) return true; + if (!pmVisualInitWindow (&kapa, "PPSub:Images")) { + isVisual = false; + return false; + } + + pmVisualScaleImage(kapa, image, "Image", 0, true); + pmVisualScaleImage(kapa, ref, "Reference", 1, true); + pmVisualScaleImage(kapa, sub, "Subtraction", 2, true); + pmVisualAskUser(&plotImage, &isVisual); + return true; +} + +static bool plotStampLocations(pmSubtractionStampList *stamps, pmReadout *ro) { + + if (!pmVisualScaleImage(kapa2, ro->image, "Stamp_master_image", 0, true)) { + fprintf(stderr, "Cannot display postage stamp master image. Skipping \n"); + return false; + } + + int Noverlay; + KiiOverlay *overlay; + + // note: this uses the Ohana allocation tools: + // ALLOCATE (overlay, KiiOverlay, 3*peaks->n + 1); + ALLOCATE (overlay, KiiOverlay, stamps->num); + + Noverlay = 0; + char boxID[10]; + for (int i = 0; i < stamps->num; i++) { + + pmSubtractionStamp *stamp = stamps->stamps->data[i]; + if (stamp == NULL) continue; + + overlay[Noverlay].type = KII_OVERLAY_BOX; + overlay[Noverlay].x = stamp->x; + overlay[Noverlay].y = stamp->y; + overlay[Noverlay].dx = 40.0; + overlay[Noverlay].dy = 40.0; + overlay[Noverlay].angle = 0.0; + sprintf(boxID, "%d", i); + overlay[Noverlay].text = boxID; + Noverlay ++; + } + + KiiLoadOverlay (kapa2, overlay, Noverlay, "red"); + FREE (overlay); + return true; +} + +#else +bool pmSubtractionSetVisual (bool mode) {return true;} +bool pmSubtractionVisualClose() {return true;} +bool pmSubtractionVisualPlotConvKernels(psImage *convKernels) {return true;} +bool pmSubtractionVisualPlotStamps(pmSubtractionStampList *stamps, pmReadout *ro) {return true;} +bool pmSubtractionVisualPlotLeastSquares(pmSubtractionStampList *stamps) {return true;} +bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub) {return true;} +#endif Index: trunk/psModules/src/imcombine/pmSubtractionVisual.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionVisual.h (revision 21422) +++ trunk/psModules/src/imcombine/pmSubtractionVisual.h (revision 21422) @@ -0,0 +1,11 @@ +#ifndef PM_SUBTRACTION_VISUAL_H +#define PM_SUBTRACTION_VISUAL_H + +bool pmSubtractionSetVisual (bool mode); +bool pmSubtractionVisualClose(); +bool pmSubtractionVisualPlotConvKernels(psImage *convKernels); +bool pmSubtractionVisualPlotStamps(pmSubtractionStampList *stamps, pmReadout *ro); +bool pmSubtractionVisualPlotLeastSquares(pmSubtractionStampList *stamps); +bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub); + +#endif Index: trunk/psModules/src/psmodules.h =================================================================== --- trunk/psModules/src/psmodules.h (revision 21387) +++ trunk/psModules/src/psmodules.h (revision 21422) @@ -9,4 +9,5 @@ #include #include +#include // XXX the following headers define constructs needed by the elements below @@ -98,4 +99,5 @@ #include #include +#include // the following headers are from psModule:objects