Index: /trunk/psModules/src/astrom/pmAstrometryObjects.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 21421) +++ /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 21421) +++ /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 21421) +++ /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 21421) +++ /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 21421) +++ /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 21421) +++ /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 21421) +++ /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 21421) +++ /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 21421) +++ /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 21421) +++ /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 Index: /trunk/psastro/src/Makefile.am =================================================================== --- /trunk/psastro/src/Makefile.am (revision 21421) +++ /trunk/psastro/src/Makefile.am (revision 21422) @@ -53,5 +53,4 @@ psastroErrorCodes.c \ psastroVersion.c \ - psastroVisual.c \ psastroDefineFiles.c \ psastroAnalysis.c \ Index: /trunk/psastro/src/psastro.h =================================================================== --- /trunk/psastro/src/psastro.h (revision 21421) +++ /trunk/psastro/src/psastro.h (revision 21422) @@ -1,8 +1,8 @@ /** @file psastro.h * - * @brief This file defines the library functions available to external - * programs. + * @brief This file defines the library functions available to external + * programs. * - * It must be included by programs which are compiled against + * It must be included by programs which are compiled against * psphot functions. * @@ -10,6 +10,6 @@ * * @author IfA - * @version $Revision: 1.48 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-07 02:03:34 $ + * @version $Revision: 1.49 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-09 21:25:34 $ * Copyright 2009 Institute for Astronomy, University of Hawaii */ @@ -30,4 +30,5 @@ # define SIGN(X) (((X) == 0) ? 0 : ((fabs((double)(X))) / (X))) +#if 0 /** * this structure represents a fit to the logN / logS curve for a set of stars @@ -43,4 +44,5 @@ int sPeak; ///< sum of stars to peak bin } pmLumFunc; +#endif bool psastroDataSave (pmConfig *config); @@ -95,18 +97,4 @@ psString psastroVersionLong(void); -// psastroVisual functions -bool psastroSetVisual (bool mode); -bool psastroVisualClose(); -bool psastroVisualPlotLuminosityFunction (psVector *lnMag, psVector *Mag, pmLumFunc *lumFunc, pmLumFunc *rawFunc); -bool psastroVisualPlotRawStars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe); -bool psastroVisualPlotRefStars (psArray *refstars, psMetadata *recipe); -bool psastroVisualPlotRemoveClumps (psArray *input, psImage *count, int scale, float limit); -bool psastroVisualPlotFixChips (pmFPAfile *input, psVector *xOld, psVector *yOld); -bool psastroVisualPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe); -bool psastroVisualPlotAstromGuessCheck (psVector *cornerPo, psVector *cornerQo, psVector *cornerPn, psVector *cornerQn, psVector *cornerPd, psVector *cornerQd); -bool psastroVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe); -bool psastroVisualPlotCommonScale (pmFPA *fpa, psVector *oldScale); -bool psastroVisualPlotMosaicMatches (psArray *rawstars, psArray *refstars, psArray *match, int iteration, psMetadata *recipe); - // demo plots bool psastroPlotRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe); Index: /trunk/psastro/src/psastroArguments.c =================================================================== --- /trunk/psastro/src/psastroArguments.c (revision 21421) +++ /trunk/psastro/src/psastroArguments.c (revision 21422) @@ -6,6 +6,6 @@ * * @author IfA - * @version $Revision: 1.33 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-07 02:03:34 $ + * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-09 21:25:34 $ * Copyright 2009 Institute for Astronomy, University of Hawaii */ @@ -100,5 +100,4 @@ if ((N = psArgumentGet (argc, argv, "-visual"))) { psArgumentRemove (N, &argc, argv); - psastroSetVisual (true); pmAstromSetVisual (true); } Index: /trunk/psastro/src/psastroAstromGuess.c =================================================================== --- /trunk/psastro/src/psastroAstromGuess.c (revision 21421) +++ /trunk/psastro/src/psastroAstromGuess.c (revision 21422) @@ -6,6 +6,6 @@ * * @author IfA - * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-07 02:03:34 $ + * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-09 21:25:34 $ * Copyright 2009 Institute for Astronomy, University of Hawaii */ @@ -42,6 +42,6 @@ psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE); if (!recipe) { - psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!"); - return false; + psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!"); + return false; } @@ -49,5 +49,5 @@ bool useModel = psMetadataLookupBool (&status, config->arguments, "PSASTRO.USE.MODEL"); if (!status) { - useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL"); + useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL"); } @@ -55,6 +55,6 @@ pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); if (!input) { - psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); - return false; + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); + return false; } @@ -62,7 +62,7 @@ double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE"); if (!status) { - psError(PS_ERR_IO, true, "Failed to lookup pixel scale"); - return false; - } + psError(PS_ERR_IO, true, "Failed to lookup pixel scale"); + return false; + } psVector *cornerL = psVectorAllocEmpty (100, PS_TYPE_F32); @@ -80,5 +80,5 @@ bool bilevelAstrometry = false; if (!useModel) { - psastroAstromGuessSetFPA (fpa, &bilevelAstrometry); + psastroAstromGuessSetFPA (fpa, &bilevelAstrometry); } @@ -88,35 +88,35 @@ if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; } - if (!useModel) { - if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue; - } + if (!useModel) { + if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue; + } if (newFPA) { newFPA = false; - while (fpa->toSky->R < 0) fpa->toSky->R += 2.0*M_PI; - while (fpa->toSky->R > 2.0*M_PI) fpa->toSky->R -= 2.0*M_PI; + while (fpa->toSky->R < 0) fpa->toSky->R += 2.0*M_PI; + while (fpa->toSky->R > 2.0*M_PI) fpa->toSky->R -= 2.0*M_PI; RAminSky = fpa->toSky->R - M_PI; RAmaxSky = fpa->toSky->R + M_PI; } - // report and save the current best guess for the chip 0,0 pixel coordinates - { - psPlane ptCH, ptFP, ptTP; - psSphere ptSky; - - ptCH.x = 0; - ptCH.y = 0; - psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); - psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP); - psDeproject (&ptSky, &ptTP, fpa->toSky); - psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d); - - psVectorAppend (cornerL, ptFP.x); - psVectorAppend (cornerM, ptFP.y); - psVectorAppend (cornerP, ptTP.x); - psVectorAppend (cornerQ, ptTP.y); - psVectorAppend (cornerR, ptSky.r); - psVectorAppend (cornerD, ptSky.d); - } + // report and save the current best guess for the chip 0,0 pixel coordinates + { + psPlane ptCH, ptFP, ptTP; + psSphere ptSky; + + ptCH.x = 0; + ptCH.y = 0; + psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); + psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP); + psDeproject (&ptSky, &ptTP, fpa->toSky); + psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d); + + psVectorAppend (cornerL, ptFP.x); + psVectorAppend (cornerM, ptFP.y); + psVectorAppend (cornerP, ptTP.x); + psVectorAppend (cornerQ, ptTP.y); + psVectorAppend (cornerR, ptSky.r); + psVectorAppend (cornerD, ptSky.d); + } // apply the new WCS guess data to all of the data in the readouts @@ -132,5 +132,5 @@ if (rawstars == NULL) { continue; } - *nStars += rawstars->n; + *nStars += rawstars->n; for (int i = 0; i < rawstars->n; i++) { pmAstromObj *raw = rawstars->data[i]; @@ -151,14 +151,14 @@ } - // dump or plot the resulting projected positions - if (psTraceGetLevel("psastro.dump") > 0) { - psastroDumpRawstars (rawstars, fpa, chip); - } - - psastroVisualPlotRawStars(rawstars, fpa, chip, recipe); - - if (psTraceGetLevel("psastro.plot") > 0) { - psastroPlotRawstars (rawstars, fpa, chip, recipe); - } + // dump or plot the resulting projected positions + if (psTraceGetLevel("psastro.dump") > 0) { + psastroDumpRawstars (rawstars, fpa, chip); + } + + pmAstromVisualPlotRawStars(rawstars, fpa, chip, recipe); + + if (psTraceGetLevel("psastro.plot") > 0) { + psastroPlotRawstars (rawstars, fpa, chip, recipe); + } } } @@ -170,12 +170,12 @@ psMetadataAddS32 (recipe, PS_LIST_TAIL, "NTOTSTAR", PS_META_REPLACE, "", *nStars); if (*nStars == 0) { - psLogMsg ("psastro", 2, "no sources available for astrometry\n"); - psFree (view); - return true; - } - - psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n", - DEG_RAD*RAmin, DEG_RAD*DECmin, - DEG_RAD*RAmax, DEG_RAD*DECmax); + psLogMsg ("psastro", 2, "no sources available for astrometry\n"); + psFree (view); + return true; + } + + psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n", + DEG_RAD*RAmin, DEG_RAD*DECmin, + DEG_RAD*RAmax, DEG_RAD*DECmax); psMetadataAddF32 (recipe, PS_LIST_TAIL, "RA_MIN", PS_META_REPLACE, "", RAmin); @@ -216,13 +216,13 @@ pmHDU *hdu = pmFPAviewThisHDU (view, fpa); if (bilevelAstrometry) { - if (!pmAstromReadBilevelChip (chip, hdu->header)) { - psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); - return false; - } + if (!pmAstromReadBilevelChip (chip, hdu->header)) { + psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); + return false; + } } else { - if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) { - psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); - return false; - } + if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) { + psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); + return false; + } } return true; @@ -238,12 +238,12 @@ // load mosaic-level astrometry? if (phu) { - char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1"); - if (ctype) { - *bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); - } + char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1"); + if (ctype) { + *bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); + } } if (*bilevelAstrometry) { - pmAstromReadBilevelMosaic (fpa, phu->header); - } + pmAstromReadBilevelMosaic (fpa, phu->header); + } psFree (view); return true; @@ -258,6 +258,6 @@ pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); if (!input) { - psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); - return false; + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); + return false; } @@ -290,55 +290,55 @@ if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; } - // XXX we are currently inconsistent with marking the good vs the bad data - // psastroChipAstrom sets data_exists to false if the fit is bad. this is - // probably wrong since it implies there is no data! - - // skip chips for which the astrometry failed (NASTRO == 0) - if (!chip->cells->n) goto skip_chip; - pmCell *cell = chip->cells->data[0]; - if (!cell) goto skip_chip; - - if (!cell->readouts->n) goto skip_chip; - pmReadout *readout = cell->readouts->data[0]; - if (!readout) goto skip_chip; - - psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"); - if (!updates) goto skip_chip; - - int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO"); - if (!nAstro) goto skip_chip; - - float astError = psMetadataLookupF32 (&status, updates, "CERROR"); - if (fabs(astError) < 1e-6) goto skip_chip; - - psPlane ptCH, ptFP, ptTP; - psSphere ptSky; - - ptCH.x = 0; - ptCH.y = 0; - psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); - psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP); - psDeproject (&ptSky, &ptTP, fpa->toSky); - psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d); - - // new corner locations based on the calibrated astrometry - psVectorAppend (cornerLn, ptFP.x); - psVectorAppend (cornerMn, ptFP.y); - psVectorAppend (cornerPn, ptTP.x); - psVectorAppend (cornerQn, ptTP.y); - psVectorAppend (cornerRn, ptSky.r); - psVectorAppend (cornerDn, ptSky.d); - psVectorAppend (cornerMK, 0); - continue; + // XXX we are currently inconsistent with marking the good vs the bad data + // psastroChipAstrom sets data_exists to false if the fit is bad. this is + // probably wrong since it implies there is no data! + + // skip chips for which the astrometry failed (NASTRO == 0) + if (!chip->cells->n) goto skip_chip; + pmCell *cell = chip->cells->data[0]; + if (!cell) goto skip_chip; + + if (!cell->readouts->n) goto skip_chip; + pmReadout *readout = cell->readouts->data[0]; + if (!readout) goto skip_chip; + + psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"); + if (!updates) goto skip_chip; + + int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO"); + if (!nAstro) goto skip_chip; + + float astError = psMetadataLookupF32 (&status, updates, "CERROR"); + if (fabs(astError) < 1e-6) goto skip_chip; + + psPlane ptCH, ptFP, ptTP; + psSphere ptSky; + + ptCH.x = 0; + ptCH.y = 0; + psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); + psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP); + psDeproject (&ptSky, &ptTP, fpa->toSky); + psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d); + + // new corner locations based on the calibrated astrometry + psVectorAppend (cornerLn, ptFP.x); + psVectorAppend (cornerMn, ptFP.y); + psVectorAppend (cornerPn, ptTP.x); + psVectorAppend (cornerQn, ptTP.y); + psVectorAppend (cornerRn, ptSky.r); + psVectorAppend (cornerDn, ptSky.d); + psVectorAppend (cornerMK, 0); + continue; skip_chip: - // new corner locations based on the calibrated astrometry - psVectorAppend (cornerLn, 0.0); - psVectorAppend (cornerMn, 0.0); - psVectorAppend (cornerPn, 0.0); - psVectorAppend (cornerQn, 0.0); - psVectorAppend (cornerRn, 0.0); - psVectorAppend (cornerDn, 0.0); - psVectorAppend (cornerMK, 1); + // new corner locations based on the calibrated astrometry + psVectorAppend (cornerLn, 0.0); + psVectorAppend (cornerMn, 0.0); + psVectorAppend (cornerPn, 0.0); + psVectorAppend (cornerQn, 0.0); + psVectorAppend (cornerRn, 0.0); + psVectorAppend (cornerDn, 0.0); + psVectorAppend (cornerMK, 1); } @@ -349,14 +349,14 @@ for (int i = 0; i < cornerRo->n; i++) { - - psPlane ptTP; - psSphere ptSky; - - ptSky.r = cornerRo->data.F32[i]; - ptSky.d = cornerDo->data.F32[i]; - - psProject (&ptTP, &ptSky, fpa->toSky); - psVectorAppend (cornerPs, ptTP.x); - psVectorAppend (cornerQs, ptTP.y); + + psPlane ptTP; + psSphere ptSky; + + ptSky.r = cornerRo->data.F32[i]; + ptSky.d = cornerDo->data.F32[i]; + + psProject (&ptTP, &ptSky, fpa->toSky); + psVectorAppend (cornerPs, ptTP.x); + psVectorAppend (cornerQs, ptTP.y); } @@ -364,9 +364,9 @@ map->x->coeffMask[1][1] = PS_POLY_MASK_SET; map->y->coeffMask[1][1] = PS_POLY_MASK_SET; - + // fit the valid chips, mask the invalid chips psVectorFitPolynomial2D (map->x, cornerMK, 1, cornerPn, NULL, cornerPs, cornerQs); psVectorFitPolynomial2D (map->y, cornerMK, 1, cornerQn, NULL, cornerPs, cornerQs); - + // apply the linear fit... psVector *cornerPf = psPolynomial2DEvalVector (map->x, cornerPs, cornerQs); @@ -377,5 +377,5 @@ psVector *cornerQd = (psVector *) psBinaryOp (NULL, cornerQn, "-", cornerQf); - psastroVisualPlotAstromGuessCheck (cornerPo, cornerQo, cornerPn, cornerQn, cornerPd, cornerQd); + pmAstromVisualPlotAstromGuessCheck (cornerPo, cornerQo, cornerPn, cornerQn, cornerPd, cornerQd); psStats *statsP = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); @@ -387,5 +387,5 @@ float angle = atan2 (map->y->coeff[1][0], map->x->coeff[1][0]); float scale = hypot (map->y->coeff[1][0], map->x->coeff[1][0]); - + psLogMsg ("psastro", 3, "boresite offset : %f,%f\n", map->x->coeff[0][0], map->y->coeff[0][0]); psLogMsg ("psastro", 3, "boresite angle : %f, scale: %f", angle*PS_DEG_RAD, scale); @@ -395,7 +395,7 @@ psMetadata *header = psMetadataLookupMetadata (&status, input->fpa->analysis, "PSASTRO.HEADER"); if (!header) { - header = psMetadataAlloc(); - psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_META_REPLACE, "psastro header stats", header); - psFree (header); // drop this reference + header = psMetadataAlloc(); + psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_META_REPLACE, "psastro header stats", header); + psFree (header); // drop this reference } @@ -408,11 +408,11 @@ if (DEBUG) { - FILE *f = fopen ("corners.dat", "w"); - for (int i = 0; i < cornerRo->n; i++) { - fprintf (f, "%10.6f %10.6f %9.2f %9.2f %9.2f %9.2f | %10.6f %10.6f %9.2f %9.2f %9.2f %9.2f\n", - cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i], - cornerRo->data.F32[i], cornerDo->data.F32[i], cornerPo->data.F32[i], cornerQo->data.F32[i], cornerLo->data.F32[i], cornerMo->data.F32[i]); - } - fclose (f); + FILE *f = fopen ("corners.dat", "w"); + for (int i = 0; i < cornerRo->n; i++) { + fprintf (f, "%10.6f %10.6f %9.2f %9.2f %9.2f %9.2f | %10.6f %10.6f %9.2f %9.2f %9.2f %9.2f\n", + cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i], + cornerRo->data.F32[i], cornerDo->data.F32[i], cornerPo->data.F32[i], cornerQo->data.F32[i], cornerLo->data.F32[i], cornerMo->data.F32[i]); + } + fclose (f); } @@ -437,5 +437,5 @@ psFree (map); psFree (view); - + return true; Index: /trunk/psastro/src/psastroCleanup.c =================================================================== --- /trunk/psastro/src/psastroCleanup.c (revision 21421) +++ /trunk/psastro/src/psastroCleanup.c (revision 21422) @@ -6,6 +6,6 @@ * * @author IfA - * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-07 02:03:34 $ + * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-09 21:25:34 $ * Copyright 2009 Institute for Astronomy, University of Hawaii */ @@ -16,4 +16,5 @@ psFree (config); + pmAstromVisualClose (); psTimerStop (); @@ -23,6 +24,4 @@ pmConceptsDone (); pmConfigDone (); - psastroVisualClose (); - pmAstromVisualClose (); fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psastro"); // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psastro"); Index: /trunk/psastro/src/psastroFixChips.c =================================================================== --- /trunk/psastro/src/psastroFixChips.c (revision 21421) +++ /trunk/psastro/src/psastroFixChips.c (revision 21422) @@ -1,11 +1,11 @@ /** @file psastroFixChips.c * - * @brief + * @brief * * @ingroup libpsastro * * @author IfA - * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-07 02:03:34 $ + * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-09 21:25:34 $ * Copyright 2009 Institute for Astronomy, University of Hawaii */ @@ -25,5 +25,5 @@ bool fixChips = psMetadataLookupBool (&status, config->arguments, "PSASTRO.FIX.CHIPS"); if (!status) { - fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS"); + fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS"); } if (!fixChips) return true; @@ -37,6 +37,6 @@ pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); if (!input) { - psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); - return false; + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); + return false; } @@ -60,6 +60,6 @@ // files associated with the science image if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { - psError (PS_ERR_IO, false, "Can't load the astrometry model file"); - return false; + psError (PS_ERR_IO, false, "Can't load the astrometry model file"); + return false; } @@ -76,100 +76,100 @@ if (DEBUG) { - f = fopen ("corners.raw.dat", "w"); - chipName = NULL; + f = fopen ("corners.raw.dat", "w"); + chipName = NULL; } pmChip *obsChip = NULL; while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { - if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; } - - // XXX we are currently inconsistent with marking the good vs the bad data - // psastroChipAstrom sets data_exists to false if the fit is bad. this is - // probably wrong since it implies there is no data! - - // skip chips for which the astrometry failed (NASTRO == 0) - if (!obsChip->cells->n) continue; - pmCell *cell = obsChip->cells->data[0]; - if (!cell) continue; - - if (!cell->readouts->n) continue; - pmReadout *readout = cell->readouts->data[0]; - if (!readout) continue; - - psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"); - if (!updates) continue; - - int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO"); - if (!nAstro) continue; - - // set the chip astrometry using the astrom file - pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa); - - psRegion *region = pmChipPixels (obsChip); - psPlane ptCP, ptFP; - - ptCP.x = region->x0; ptCP.y = region->y0; - psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); - xObs->data.F32[nPts] = ptFP.x; - yObs->data.F32[nPts] = ptFP.y; - psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP); - xRef->data.F32[nPts] = ptFP.x; - yRef->data.F32[nPts] = ptFP.y; - - if (DEBUG) { - chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME"); - fprintf (f, "%s %f %f %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]); - } - nPts ++; - - ptCP.x = region->x0; ptCP.y = region->y1; - psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); - xObs->data.F32[nPts] = ptFP.x; - yObs->data.F32[nPts] = ptFP.y; - psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP); - xRef->data.F32[nPts] = ptFP.x; - yRef->data.F32[nPts] = ptFP.y; - - if (DEBUG) { - chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME"); - fprintf (f, "%s %f %f %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]); - } - nPts ++; - - ptCP.x = region->x1; ptCP.y = region->y1; - psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); - xObs->data.F32[nPts] = ptFP.x; - yObs->data.F32[nPts] = ptFP.y; - psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP); - xRef->data.F32[nPts] = ptFP.x; - yRef->data.F32[nPts] = ptFP.y; - - if (DEBUG) { - chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME"); - fprintf (f, "%s %f %f %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]); - } - nPts ++; - - ptCP.x = region->x1; ptCP.y = region->y0; - psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); - xObs->data.F32[nPts] = ptFP.x; - yObs->data.F32[nPts] = ptFP.y; - psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP); - xRef->data.F32[nPts] = ptFP.x; - yRef->data.F32[nPts] = ptFP.y; - - if (DEBUG) { - chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME"); - fprintf (f, "%s %f %f %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]); - } - nPts ++; - - psFree (region); + if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; } + + // XXX we are currently inconsistent with marking the good vs the bad data + // psastroChipAstrom sets data_exists to false if the fit is bad. this is + // probably wrong since it implies there is no data! + + // skip chips for which the astrometry failed (NASTRO == 0) + if (!obsChip->cells->n) continue; + pmCell *cell = obsChip->cells->data[0]; + if (!cell) continue; + + if (!cell->readouts->n) continue; + pmReadout *readout = cell->readouts->data[0]; + if (!readout) continue; + + psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"); + if (!updates) continue; + + int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO"); + if (!nAstro) continue; + + // set the chip astrometry using the astrom file + pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa); + + psRegion *region = pmChipPixels (obsChip); + psPlane ptCP, ptFP; + + ptCP.x = region->x0; ptCP.y = region->y0; + psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); + xObs->data.F32[nPts] = ptFP.x; + yObs->data.F32[nPts] = ptFP.y; + psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP); + xRef->data.F32[nPts] = ptFP.x; + yRef->data.F32[nPts] = ptFP.y; + + if (DEBUG) { + chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME"); + fprintf (f, "%s %f %f %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]); + } + nPts ++; + + ptCP.x = region->x0; ptCP.y = region->y1; + psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); + xObs->data.F32[nPts] = ptFP.x; + yObs->data.F32[nPts] = ptFP.y; + psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP); + xRef->data.F32[nPts] = ptFP.x; + yRef->data.F32[nPts] = ptFP.y; + + if (DEBUG) { + chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME"); + fprintf (f, "%s %f %f %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]); + } + nPts ++; + + ptCP.x = region->x1; ptCP.y = region->y1; + psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); + xObs->data.F32[nPts] = ptFP.x; + yObs->data.F32[nPts] = ptFP.y; + psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP); + xRef->data.F32[nPts] = ptFP.x; + yRef->data.F32[nPts] = ptFP.y; + + if (DEBUG) { + chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME"); + fprintf (f, "%s %f %f %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]); + } + nPts ++; + + ptCP.x = region->x1; ptCP.y = region->y0; + psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); + xObs->data.F32[nPts] = ptFP.x; + yObs->data.F32[nPts] = ptFP.y; + psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP); + xRef->data.F32[nPts] = ptFP.x; + yRef->data.F32[nPts] = ptFP.y; + + if (DEBUG) { + chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME"); + fprintf (f, "%s %f %f %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]); + } + nPts ++; + + psFree (region); } xObs->n = yObs->n = xRef->n = yRef->n = nPts; if (DEBUG) fclose (f); - + psPlaneTransform *map = psPlaneTransformAlloc (1, 1); - + psVector *mask = psVectorAlloc (nPts, PS_TYPE_VECTOR_MASK); psVectorInit (mask, 0); @@ -179,25 +179,25 @@ for (int i = 0; i < 3; i++) { - psVectorClipFitPolynomial2D (map->x, stats, mask, 0xff, xObs, NULL, xRef, yRef); - psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, xObs->n); - - psVectorClipFitPolynomial2D (map->y, stats, mask, 0xff, yObs, NULL, xRef, yRef); - psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, yObs->n); + psVectorClipFitPolynomial2D (map->x, stats, mask, 0xff, xObs, NULL, xRef, yRef); + psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, xObs->n); + + psVectorClipFitPolynomial2D (map->y, stats, mask, 0xff, yObs, NULL, xRef, yRef); + psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, yObs->n); } // loop over all chips, select the outliers, and replace the measured astrometry with the model - // the measured transformation above must be applied to make the comparison, and also then applied to the + // the measured transformation above must be applied to make the comparison, and also then applied to the // model transformation if (DEBUG) { - f = fopen ("corners.fit.dat", "w"); - for (int i = 0; i < xObs->n; i++) { - psPlane obsCoord, refCoord; - refCoord.x = xRef->data.F32[i]; - refCoord.y = yRef->data.F32[i]; - psPlaneTransformApply (&obsCoord, map, &refCoord); - fprintf (f, "%f %f %f %f %f %f\n", xObs->data.F32[i], yObs->data.F32[i], xRef->data.F32[i], yRef->data.F32[i], obsCoord.x, obsCoord.y); - } - fclose (f); + f = fopen ("corners.fit.dat", "w"); + for (int i = 0; i < xObs->n; i++) { + psPlane obsCoord, refCoord; + refCoord.x = xRef->data.F32[i]; + refCoord.y = yRef->data.F32[i]; + psPlaneTransformApply (&obsCoord, map, &refCoord); + fprintf (f, "%f %f %f %f %f %f\n", xObs->data.F32[i], yObs->data.F32[i], xRef->data.F32[i], yRef->data.F32[i], obsCoord.x, obsCoord.y); + } + fclose (f); } @@ -211,125 +211,125 @@ while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { - psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process); - if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; } - - // set the chip astrometry using the astrom file - pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa); - - // bad Astrometry test: ref pixel or angle outside nominal - - psPlane refPixel = {0.0, 0.0, 0.0, 0.0}; - psPlane obsCoord, refCoord, tmpCoord; - - // find location of 0,0 pixel in focal plane coords for this chip - psPlaneTransformApply (&obsCoord, obsChip->toFPA, &refPixel); - - // find location of 0,0 pixel in focal plane coords for ref chip - // apply the global field rotation and offset before comparing - psPlaneTransformApply (&tmpCoord, refChip->toFPA, &refPixel); - psPlaneTransformApply (&refCoord, map, &tmpCoord); - - psPlane offPixel = {100.0, 0.0, 100.0, 0.0}; - psPlane obsOffPt, refOffPt; - - // find location of 0,0 pixel in focal plane coords for this chip - psPlaneTransformApply (&obsOffPt, obsChip->toFPA, &offPixel); - - // find location of 0,0 pixel in focal plane coords for ref chip - psPlaneTransformApply (&tmpCoord, refChip->toFPA, &offPixel); - psPlaneTransformApply (&refOffPt, map, &tmpCoord); - - double obsAngle = PM_DEG_RAD*atan2 (obsOffPt.y - obsCoord.y, obsOffPt.x - obsCoord.x); - double refAngle = PM_DEG_RAD*atan2 (refOffPt.y - refCoord.y, refOffPt.x - refCoord.x); - - bool badAstrom = false; - badAstrom |= fabs(obsCoord.x - refCoord.x) > pixelTol; - badAstrom |= fabs(obsCoord.y - refCoord.y) > pixelTol; - badAstrom |= fabs(obsAngle - refAngle) > angleTol; - - fprintf (stderr, "chip %d, angle: %f, pixel: %f,%f\n", view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y); - - // XXX for now, just use first readout - pmCell *cell = obsChip->cells->data[0]; - pmReadout *readout = cell->readouts->data[0]; - - // update the header (pull or create local view to entry on readout->analysis) - psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"); - if (!updates) { - updates = psMetadataAlloc (); - psMetadataAddMetadata (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_META_REPLACE, "psastro header stats", updates); - psFree (updates); - } - - psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DX", PS_META_REPLACE, "chip x offset wrt model", obsCoord.x - refCoord.x); - psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DY", PS_META_REPLACE, "chip y offset wrt model", obsCoord.y - refCoord.y); - psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DT", PS_META_REPLACE, "chip rot offset wrt model", obsAngle - refAngle); - - // for successful chips, save the measured offsets in the header - if (!badAstrom) continue; - - // XXX for now, let's just fail on the bad chips. In the future, let's try to recover, but we still need to - // catch the failures relative to the model - psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", 0); - continue; - - psLogMsg ("psastro", PS_LOG_INFO, "fixing chip %d, angle: %f, pixel: %f,%f\n", - view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y); - - psFree (obsChip->toFPA); - psFree (obsChip->fromFPA); - - // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter - // XXX this only works if toTPA is at most a linear transformation - psPlaneTransform *toFPA = psPlaneTransformAlloc(refChip->toFPA->x->nX, refChip->toFPA->x->nY); - for (int i = 0; i <= refChip->toFPA->x->nX; i++) { - for (int j = 0; j <= refChip->toFPA->x->nY; j++) { - double f1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->x->coeff[1][0]*refChip->toFPA->x->coeff[i][j]; - double f2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->x->coeff[0][1]*refChip->toFPA->y->coeff[i][j]; - toFPA->x->coeff[i][j] = f1 + f2; - - double g1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->y->coeff[1][0]*refChip->toFPA->x->coeff[i][j]; - double g2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->y->coeff[0][1]*refChip->toFPA->y->coeff[i][j]; - toFPA->y->coeff[i][j] = g1 + g2; - } - } - toFPA->x->coeff[0][0] += map->x->coeff[0][0]; - toFPA->y->coeff[0][0] += map->y->coeff[0][0]; - - psRegion *region = pmChipPixels (obsChip); - obsChip->toFPA = toFPA; - obsChip->fromFPA = psPlaneTransformInvert(NULL, obsChip->toFPA, *region, 50); - psFree (region); - - // use the new position to re-try the match fit - // select the raw objects for this readout - psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); - if (rawstars == NULL) { continue; } - - // select the raw objects for this readout - psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); - if (refstars == NULL) { continue; } - - // the absolute minimum number of stars is 4 (for order = 1) - if ((rawstars->n < 4) || (refstars->n < 4)) { - readout->data_exists = false; - psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)", - rawstars->n, refstars->n); - continue; - } - - psastroUpdateChipToFPA (input->fpa, obsChip, rawstars, refstars); - - // XXX update the header with info to reflect the failure - if (!psastroOneChipFit (input->fpa, obsChip, refstars, rawstars, recipe, updates)) { - readout->data_exists = false; - psLogMsg ("psastro", 3, "failed to find a solution\n"); - continue; - } - - pmAstromWriteWCS (updates, input->fpa, obsChip, NONLIN_TOL); - } - - psastroVisualPlotFixChips (input, xObs, yObs); + psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process); + if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; } + + // set the chip astrometry using the astrom file + pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa); + + // bad Astrometry test: ref pixel or angle outside nominal + + psPlane refPixel = {0.0, 0.0, 0.0, 0.0}; + psPlane obsCoord, refCoord, tmpCoord; + + // find location of 0,0 pixel in focal plane coords for this chip + psPlaneTransformApply (&obsCoord, obsChip->toFPA, &refPixel); + + // find location of 0,0 pixel in focal plane coords for ref chip + // apply the global field rotation and offset before comparing + psPlaneTransformApply (&tmpCoord, refChip->toFPA, &refPixel); + psPlaneTransformApply (&refCoord, map, &tmpCoord); + + psPlane offPixel = {100.0, 0.0, 100.0, 0.0}; + psPlane obsOffPt, refOffPt; + + // find location of 0,0 pixel in focal plane coords for this chip + psPlaneTransformApply (&obsOffPt, obsChip->toFPA, &offPixel); + + // find location of 0,0 pixel in focal plane coords for ref chip + psPlaneTransformApply (&tmpCoord, refChip->toFPA, &offPixel); + psPlaneTransformApply (&refOffPt, map, &tmpCoord); + + double obsAngle = PM_DEG_RAD*atan2 (obsOffPt.y - obsCoord.y, obsOffPt.x - obsCoord.x); + double refAngle = PM_DEG_RAD*atan2 (refOffPt.y - refCoord.y, refOffPt.x - refCoord.x); + + bool badAstrom = false; + badAstrom |= fabs(obsCoord.x - refCoord.x) > pixelTol; + badAstrom |= fabs(obsCoord.y - refCoord.y) > pixelTol; + badAstrom |= fabs(obsAngle - refAngle) > angleTol; + + fprintf (stderr, "chip %d, angle: %f, pixel: %f,%f\n", view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y); + + // XXX for now, just use first readout + pmCell *cell = obsChip->cells->data[0]; + pmReadout *readout = cell->readouts->data[0]; + + // update the header (pull or create local view to entry on readout->analysis) + psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"); + if (!updates) { + updates = psMetadataAlloc (); + psMetadataAddMetadata (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_META_REPLACE, "psastro header stats", updates); + psFree (updates); + } + + psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DX", PS_META_REPLACE, "chip x offset wrt model", obsCoord.x - refCoord.x); + psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DY", PS_META_REPLACE, "chip y offset wrt model", obsCoord.y - refCoord.y); + psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DT", PS_META_REPLACE, "chip rot offset wrt model", obsAngle - refAngle); + + // for successful chips, save the measured offsets in the header + if (!badAstrom) continue; + + // XXX for now, let's just fail on the bad chips. In the future, let's try to recover, but we still need to + // catch the failures relative to the model + psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", 0); + continue; + + psLogMsg ("psastro", PS_LOG_INFO, "fixing chip %d, angle: %f, pixel: %f,%f\n", + view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y); + + psFree (obsChip->toFPA); + psFree (obsChip->fromFPA); + + // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter + // XXX this only works if toTPA is at most a linear transformation + psPlaneTransform *toFPA = psPlaneTransformAlloc(refChip->toFPA->x->nX, refChip->toFPA->x->nY); + for (int i = 0; i <= refChip->toFPA->x->nX; i++) { + for (int j = 0; j <= refChip->toFPA->x->nY; j++) { + double f1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->x->coeff[1][0]*refChip->toFPA->x->coeff[i][j]; + double f2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->x->coeff[0][1]*refChip->toFPA->y->coeff[i][j]; + toFPA->x->coeff[i][j] = f1 + f2; + + double g1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->y->coeff[1][0]*refChip->toFPA->x->coeff[i][j]; + double g2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->y->coeff[0][1]*refChip->toFPA->y->coeff[i][j]; + toFPA->y->coeff[i][j] = g1 + g2; + } + } + toFPA->x->coeff[0][0] += map->x->coeff[0][0]; + toFPA->y->coeff[0][0] += map->y->coeff[0][0]; + + psRegion *region = pmChipPixels (obsChip); + obsChip->toFPA = toFPA; + obsChip->fromFPA = psPlaneTransformInvert(NULL, obsChip->toFPA, *region, 50); + psFree (region); + + // use the new position to re-try the match fit + // select the raw objects for this readout + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); + if (rawstars == NULL) { continue; } + + // select the raw objects for this readout + psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); + if (refstars == NULL) { continue; } + + // the absolute minimum number of stars is 4 (for order = 1) + if ((rawstars->n < 4) || (refstars->n < 4)) { + readout->data_exists = false; + psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)", + rawstars->n, refstars->n); + continue; + } + + psastroUpdateChipToFPA (input->fpa, obsChip, rawstars, refstars); + + // XXX update the header with info to reflect the failure + if (!psastroOneChipFit (input->fpa, obsChip, refstars, rawstars, recipe, updates)) { + readout->data_exists = false; + psLogMsg ("psastro", 3, "failed to find a solution\n"); + continue; + } + + pmAstromWriteWCS (updates, input->fpa, obsChip, NONLIN_TOL); + } + + pmAstromVisualPlotFixChips (input, xObs, yObs); psFree (xObs); psFree (yObs); Index: /trunk/psastro/src/psastroLoadRefstars.c =================================================================== --- /trunk/psastro/src/psastroLoadRefstars.c (revision 21421) +++ /trunk/psastro/src/psastroLoadRefstars.c (revision 21422) @@ -6,6 +6,6 @@ * * @author IfA - * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-07 02:03:34 $ + * @version $Revision: 1.36 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-09 21:25:34 $ * Copyright 2009 Institute for Astronomy, University of Hawaii */ @@ -150,5 +150,5 @@ } - psastroVisualPlotRefStars (refstars, recipe); + pmAstromVisualPlotRefStars (refstars, recipe); if (psTraceGetLevel("psastro.plot") > 0) { @@ -251,7 +251,7 @@ pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); if (!input) { - psLogMsg ("psastro", PS_LOG_DETAIL, "no supplied reference header data"); - photcode = psStringCopy ("NONE"); - return photcode; + psLogMsg ("psastro", PS_LOG_DETAIL, "no supplied reference header data"); + photcode = psStringCopy ("NONE"); + return photcode; } assert (input->fpa); @@ -259,6 +259,6 @@ *maxRho = psMetadataLookupF32(&status, recipe, "DVO.GETSTAR.MAX.RHO"); if (!status) { - psError(PSASTRO_ERR_CONFIG, false, "DVO.GETSTAR.MAX.RHO missing from recipe"); - return NULL; + psError(PSASTRO_ERR_CONFIG, false, "DVO.GETSTAR.MAX.RHO missing from recipe"); + return NULL; } @@ -278,39 +278,39 @@ if (!status) ESCAPE ("missing DVO.GETSTAR.MIN.MAG.INST"); - // PHOTCODE.DATA is a multi of metadata items + // PHOTCODE.DATA is a multi of metadata items psListIterator *iter = psListIteratorAlloc(item->data.list, PS_LIST_HEAD, false); psMetadataItem *refItem = NULL; while ((refItem = psListGetAndIncrement (iter))) { - if (refItem->type != PS_DATA_METADATA) ESCAPE ("PHOTCODE.DATA entry is not a metadata folder"); - - char *refFilter = psMetadataLookupStr (&status, refItem->data.md, "FILTER"); - if (!status) { - psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER"); - continue; - } - - // does this entry match the current filter? - if (strcmp (refFilter, filter)) continue; - - psLogMsg ("psastro", PS_LOG_DETAIL, "PHOTCODE.DATA found for filter %s", filter); - - float zeropt = psMetadataLookupF32 (&status, refItem->data.md, "ZEROPT"); - if (!status) { - psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER"); - continue; - } - photcode = psMetadataLookupStr (&status, refItem->data.md, "PHOTCODE"); - if (!status) { - psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER"); - continue; - } - - // convert the minInst to a calibrated minimum magnitude - *minMag = minInst + 2.5*log10(exptime) + zeropt; - - psFree (iter); - return photcode; - } + if (refItem->type != PS_DATA_METADATA) ESCAPE ("PHOTCODE.DATA entry is not a metadata folder"); + + char *refFilter = psMetadataLookupStr (&status, refItem->data.md, "FILTER"); + if (!status) { + psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER"); + continue; + } + + // does this entry match the current filter? + if (strcmp (refFilter, filter)) continue; + + psLogMsg ("psastro", PS_LOG_DETAIL, "PHOTCODE.DATA found for filter %s", filter); + + float zeropt = psMetadataLookupF32 (&status, refItem->data.md, "ZEROPT"); + if (!status) { + psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER"); + continue; + } + photcode = psMetadataLookupStr (&status, refItem->data.md, "PHOTCODE"); + if (!status) { + psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER"); + continue; + } + + // convert the minInst to a calibrated minimum magnitude + *minMag = minInst + 2.5*log10(exptime) + zeropt; + + psFree (iter); + return photcode; + } psFree (iter); @@ -318,5 +318,5 @@ photcode = psMetadataLookupStr(NULL, recipe, "DVO.GETSTAR.PHOTCODE"); PS_ASSERT (photcode, NULL); - + // give up and use fixed value *minMag = psMetadataLookupF32(NULL, recipe, "DVO.GETSTAR.MIN.MAG"); Index: /trunk/psastro/src/psastroLuminosityFunction.c =================================================================== --- /trunk/psastro/src/psastroLuminosityFunction.c (revision 21421) +++ /trunk/psastro/src/psastroLuminosityFunction.c (revision 21422) @@ -6,6 +6,6 @@ * * @author IfA - * @version $Revision: 1.18 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-07 02:03:34 $ + * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-09 21:25:34 $ * Copyright 2009 Institute for Astronomy, University of Hawaii */ @@ -141,5 +141,5 @@ lumFunc->sPeak = sPeak; - psastroVisualPlotLuminosityFunction(lnMag, Mag, lumFunc, rawFunc); + pmAstromVisualPlotLuminosityFunction(lnMag, Mag, lumFunc, rawFunc); psFree (lnMag); Index: /trunk/psastro/src/psastroMosaicOneChip.c =================================================================== --- /trunk/psastro/src/psastroMosaicOneChip.c (revision 21421) +++ /trunk/psastro/src/psastroMosaicOneChip.c (revision 21422) @@ -6,6 +6,6 @@ * * @author IfA - * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-07 02:03:34 $ + * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-09 21:25:34 $ * Copyright 2009 Institute for Astronomy, University of Hawaii */ @@ -160,5 +160,5 @@ //plot results - psastroVisualPlotMosaicOneChip(rawstars, refstars, match, recipe); + pmAstromVisualPlotMosaicOneChip(rawstars, refstars, match, recipe); psFree (fitStats); Index: /trunk/psastro/src/psastroMosaicSetMatch.c =================================================================== --- /trunk/psastro/src/psastroMosaicSetMatch.c (revision 21421) +++ /trunk/psastro/src/psastroMosaicSetMatch.c (revision 21422) @@ -6,6 +6,6 @@ * * @author IfA - * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-07 02:03:34 $ + * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-09 21:25:34 $ * Copyright 2009 Institute for Astronomy, University of Hawaii */ @@ -69,5 +69,5 @@ psTrace ("psastro", 4, "Matched %ld refstars\n", matches->n); - psastroVisualPlotMosaicMatches(rawstars, refstars, matches, iteration, recipe); + pmAstromVisualPlotMosaicMatches(rawstars, refstars, matches, iteration, recipe); // XXX drop the old one Index: /trunk/psastro/src/psastroOneChipFit.c =================================================================== --- /trunk/psastro/src/psastroOneChipFit.c (revision 21421) +++ /trunk/psastro/src/psastroOneChipFit.c (revision 21422) @@ -6,6 +6,6 @@ * * @author IfA - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-07 02:03:34 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-09 21:25:34 $ * Copyright 2009 Institute for Astronomy, University of Hawaii */ @@ -23,11 +23,11 @@ // default value for match/fit : radius is in pixels - REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32); + REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32); // run the match/fit sequence NITER times - REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32); + REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32); // correct radius to FP units (physical pixel scale in microns per pixel) - REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32); + REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32); RADIUS *= pixelScale; @@ -44,79 +44,79 @@ for (int iter = 0; iter < nIter; iter++) { - - char name[128]; - sprintf (name, "PSASTRO.MATCH.RADIUS.N%d", iter); - float radius = psMetadataLookupF32 (&status, recipe, name); - radius *= pixelScale; - if (!status || (radius == 0.0)) { - radius = RADIUS; - } + char name[128]; + + sprintf (name, "PSASTRO.MATCH.RADIUS.N%d", iter); + float radius = psMetadataLookupF32 (&status, recipe, name); + radius *= pixelScale; + if (!status || (radius == 0.0)) { + radius = RADIUS; + } - // use small radius to match stars - match = pmAstromRadiusMatchFP (rawstars, refstars, radius); - if (match == NULL) { - psLogMsg ("psastro", 3, "failed to find radius-matched sources\n"); - return false; - } + // use small radius to match stars + match = pmAstromRadiusMatchFP (rawstars, refstars, radius); + if (match == NULL) { + psLogMsg ("psastro", 3, "failed to find radius-matched sources\n"); + return false; + } - // modify the order to correspond to the actual number of matched stars: - int Ndof_min = 3; - int order_max = 0.5*(3 + sqrt(4*match->n - 4*Ndof_min + 1)); - order = PS_MIN (order, order_max); + // modify the order to correspond to the actual number of matched stars: + int Ndof_min = 3; + int order_max = 0.5*(3 + sqrt(4*match->n - 4*Ndof_min + 1)); + order = PS_MIN (order, order_max); - // if ((match->n < 11) && (order >= 3)) order = 2; - // if ((match->n < 7) && (order >= 2)) order = 1; - // if ((match->n < 4) && (order >= 1)) order = 0; + // if ((match->n < 11) && (order >= 3)) order = 2; + // if ((match->n < 7) && (order >= 2)) order = 1; + // if ((match->n < 4) && (order >= 1)) order = 0; - if (order < 1) { - psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n); - psFree (match); - return false; - } + if (order < 1) { + psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n); + psFree (match); + return false; + } - // create output toFPA; set masks appropriate to the Elixir DVO astrometry format - psFree (chip->toFPA); - chip->toFPA = psPlaneTransformAlloc (order, order); - for (int i = 0; i <= chip->toFPA->x->nX; i++) { - for (int j = 0; j <= chip->toFPA->x->nY; j++) { - if (i + j > order) { - chip->toFPA->x->coeffMask[i][j] = PS_POLY_MASK_SET; - chip->toFPA->y->coeffMask[i][j] = PS_POLY_MASK_SET; - } - } - } + // create output toFPA; set masks appropriate to the Elixir DVO astrometry format + psFree (chip->toFPA); + chip->toFPA = psPlaneTransformAlloc (order, order); + for (int i = 0; i <= chip->toFPA->x->nX; i++) { + for (int j = 0; j <= chip->toFPA->x->nY; j++) { + if (i + j > order) { + chip->toFPA->x->coeffMask[i][j] = PS_POLY_MASK_SET; + chip->toFPA->y->coeffMask[i][j] = PS_POLY_MASK_SET; + } + } + } - // XXX allow statitic to be set by the user - // fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); - fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA"); - fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER"); + // XXX allow statitic to be set by the user + // fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); + fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA"); + fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER"); - // improved fit for astrometric terms - results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats); - if (!results) { - psLogMsg ("psastro", 3, "failed to perform the matched fit\n"); - psFree (match); - psFree (fitStats); - return false; - } - - // determine fromFPA transformation and apply new transformation to raw & ref stars - psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); - - // toSky converts from FPA & TPA units (microns) to sky units (radians) - float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD; - // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale; - float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale; - int astNstar = results->yStats->clippedNvalues; - psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar); + // improved fit for astrometric terms + results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats); + if (!results) { + psLogMsg ("psastro", 3, "failed to perform the matched fit\n"); + psFree (match); + psFree (fitStats); + return false; + } - if (iter < nIter - 1) { - psFree (fitStats); - psFree (results); - psFree (match); - } + // determine fromFPA transformation and apply new transformation to raw & ref stars + psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); + + // toSky converts from FPA & TPA units (microns) to sky units (radians) + float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD; + // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale; + float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale; + int astNstar = results->yStats->clippedNvalues; + psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar); + + if (iter < nIter - 1) { + psFree (fitStats); + psFree (results); + psFree (match); + } } @@ -139,9 +139,9 @@ if (astError > maxError) { psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError); - validSolution = false; + validSolution = false; } if (astNstar < minNstar) { psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar); - validSolution = false; + validSolution = false; } @@ -150,9 +150,9 @@ psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR", PS_META_REPLACE, "astrometry error (arcsec)", astError); if (validSolution) { - psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar)); - psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", astNstar); + psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar)); + psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", astNstar); } else { - psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0); - psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", 0); + psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0); + psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", 0); } psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX", PS_META_REPLACE, "equinox of ref catalog", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars @@ -160,16 +160,16 @@ // XXX drop from here : determine fromFPA transformation and apply new transformation to raw & ref stars // psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); - + // XXX check if we correctly applied the new transformation: if (psTraceGetLevel("psastro.dump") > 0) { - psastroDumpRawstars (rawstars, fpa, chip); - psastroDumpMatchedStars ("match.dat", rawstars, refstars, match); - psastroDumpStars (refstars, "refstars.cal.dat"); + psastroDumpRawstars (rawstars, fpa, chip); + psastroDumpMatchedStars ("match.dat", rawstars, refstars, match); + psastroDumpStars (refstars, "refstars.cal.dat"); } - psastroVisualPlotOneChipFit (rawstars, refstars, match, recipe); + pmAstromVisualPlotOneChipFit (rawstars, refstars, match, recipe); if (psTraceGetLevel("psastro.plot") > 0) { - psastroPlotOneChipFit (rawstars, refstars, match, recipe); + psastroPlotOneChipFit (rawstars, refstars, match, recipe); } Index: /trunk/psastro/src/psastroRemoveClumps.c =================================================================== --- /trunk/psastro/src/psastroRemoveClumps.c (revision 21421) +++ /trunk/psastro/src/psastroRemoveClumps.c (revision 21422) @@ -6,6 +6,6 @@ * * @author IfA - * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-07 02:03:34 $ + * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-09 21:25:34 $ * Copyright 2009 Institute for Astronomy, University of Hawaii */ @@ -69,5 +69,5 @@ psTrace ("psastro", 4, "skipping stars in cells with more than %f stars\n", limit); - psastroVisualPlotRemoveClumps (input, count, scale, limit); + pmAstromVisualPlotRemoveClumps (input, count, scale, limit); // find and exclude objects in bad pixels Index: /trunk/psastro/src/psastroUtils.c =================================================================== --- /trunk/psastro/src/psastroUtils.c (revision 21421) +++ /trunk/psastro/src/psastroUtils.c (revision 21422) @@ -6,6 +6,6 @@ * * @author IfA - * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-07 02:03:34 $ + * @version $Revision: 1.25 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-09 21:25:34 $ * Copyright 2009 Institute for Astronomy, University of Hawaii */ @@ -54,11 +54,11 @@ if (!chip->toFPA) { continue; } - if (chip->cells->n == 0) { continue; } - pmCell *cell = chip->cells->data[0]; + if (chip->cells->n == 0) { continue; } + pmCell *cell = chip->cells->data[0]; if (!cell->process || !cell->file_exists) { continue; } - if (cell->readouts->n == 0) { continue; } - pmReadout *readout = cell->readouts->data[0]; - if (! readout->data_exists) { continue; } + if (cell->readouts->n == 0) { continue; } + pmReadout *readout = cell->readouts->data[0]; + if (! readout->data_exists) { continue; } pixelScale1 = hypot (chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1]); @@ -100,5 +100,5 @@ psastroMosaicSetAstrom (fpa); if (!useExternal) { - psastroVisualPlotCommonScale (fpa, oldScale); + pmAstromVisualPlotCommonScale (fpa, oldScale); } psFree (oldScale);