Index: /branches/cnb_branch_20081011/ippconfig/recipes/psastro.config =================================================================== --- /branches/cnb_branch_20081011/ippconfig/recipes/psastro.config (revision 20334) +++ /branches/cnb_branch_20081011/ippconfig/recipes/psastro.config (revision 20335) @@ -75,5 +75,5 @@ PSASTRO.MAX.INST.MAG.RAW F32 0.0 # max instrumental magnitude for stars accepted for fitting -PSASTRO.MATCH.LUMFUNC BOOL FALSE +PSASTRO.MATCH.LUMFUNC BOOL TRUE # option may be MAX, MIN, or VALUE. for VALUE, look up Index: /branches/cnb_branch_20081011/psModules/src/astrom/Makefile.am =================================================================== --- /branches/cnb_branch_20081011/psModules/src/astrom/Makefile.am (revision 20334) +++ /branches/cnb_branch_20081011/psModules/src/astrom/Makefile.am (revision 20335) @@ -10,5 +10,6 @@ pmAstrometryModel.c \ pmAstrometryRefstars.c \ - pmAstrometryWCS.c + pmAstrometryWCS.c \ + pmAstrometryVisual.c pkginclude_HEADERS = \ @@ -19,5 +20,6 @@ pmAstrometryModel.h \ pmAstrometryRefstars.h \ - pmAstrometryWCS.h + pmAstrometryWCS.h \ + pmAstrometryVisual.h CLEANFILES = *~ Index: /branches/cnb_branch_20081011/psModules/src/astrom/pmAstrometryObjects.c =================================================================== --- /branches/cnb_branch_20081011/psModules/src/astrom/pmAstrometryObjects.c (revision 20334) +++ /branches/cnb_branch_20081011/psModules/src/astrom/pmAstrometryObjects.c (revision 20335) @@ -8,6 +8,6 @@ * @author EAM, IfA * -* @version $Revision: 1.41 $ $Name: not supported by cvs2svn $ -* @date $Date: 2008-10-10 02:33:35 $ +* @version $Revision: 1.41.4.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2008-10-22 23:35:13 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -29,5 +29,4 @@ #include // for unlink #include -#include // cnb: do I need an IFDEF? #include "pmHDU.h" @@ -35,4 +34,5 @@ #include "pmAstrometryObjects.h" #include "pmKapaPlots.h" +#include "pmAstrometryVisual.h" #define PM_ASTROMETRYOBJECTS_DEBUG 1 @@ -564,10 +564,4 @@ psF32 **D2 = gridD2->data.F32; -#ifdef PLOTS - // vectors to hold dX and dY - int nplot = raw->n * ref->n; - float dXplot[nplot]; - float dYplot[nplot]; -#endif // accumulate grids for focal plane (L,M) matches @@ -579,9 +573,4 @@ dY = ob1->FP->y - ob2->FP->y; -#ifdef PLOTS - dXplot[(i * ref->n) + j] = dX; - dYplot[(i * ref->n) + j] = dY; - // fprintf (f, "dX,dY: %8.2f %8.2f : %8.2f %8.2f : %8.2f %8.2f\n", dX, dY, ob1->FP->x, ob2->FP->x, ob1->FP->y, ob2->FP->y); -#endif // find bin coordinates for this delta-delta if (!AstromGridBin (&iX, &iY, dX, dY)) { @@ -618,5 +607,5 @@ } - # if 0 +# if 0 char line[16]; psFits *fits = psFitsOpen ("grid.image.fits", "w"); @@ -625,5 +614,5 @@ fprintf (stderr, "wrote grid image, press return to continue\n"); fgets (line, 15, stdin); - # endif +# endif // only check bins with at least 1/2 of max bin @@ -670,4 +659,7 @@ } + pmAstromVisualPlotGridMatch (raw, ref,gridNP, stats->offset.x, stats->offset.y, + maxOffpix, Scale, Offset); + psFree (imStats); // XXX EAM : This routine, and pmAstromGridMatch, need to handle failure cases better @@ -689,59 +681,4 @@ // fprintf (stderr, "sigma: nMatch: %d, nTest: %d, nTen: %d\n", stats->nMatch, stats->nTest, sort->data.U32[sort->n - 10]); -#ifdef PLOTS - /*** plot DXs and DYs ***/ - KapaSection section = {"s1", 0.00, 0.00, 1.0, 1.0}; - //KapaSection left = {"s2", 0.05, 0.05, 0.35, 0.8}; - //KapaSection right = {"s3", 0.55, 0.05, 0.35, 0.8}; - //KapaSection yprof = {"s4", 0.45, 0.05, 0.10, 0.9}; - //KapaSection xprof = {"s5", .05, .90, .35, .1}; - - Graphdata graphdata; - int kapa = pmKapaOpen(true); - KapaClearPlots(kapa); - KapaInitGraph(&graphdata); - KapaSetSection(kapa, §ion); - graphdata.xmin = stats->offset.x - 1.5 * maxOffpix; - graphdata.xmax = stats->offset.x + 1.5 * maxOffpix; - graphdata.ymin = stats->offset.y - 1.5 * maxOffpix; - graphdata.ymax = stats->offset.y + 1.5 * maxOffpix; - - KapaSetLimits(kapa, &graphdata); - KapaSetFont(kapa, "helvetica", 14); - KapaBox(kapa, &graphdata); - KapaSendLabel (kapa, "X offset", KAPA_LABEL_XM); - KapaSendLabel (kapa, "Y offset", KAPA_LABEL_YM); - KapaSendLabel (kapa, "pmAstromGridAngle residuals. Big Box: Serach Region. Small Box: Correlation Peak.", - KAPA_LABEL_XP); - graphdata.style = 2; - graphdata.ptype = 0; - graphdata.size = 0.4; - graphdata.color = KapaColorByName ("black"); - KapaPrepPlot(kapa, nplot, &graphdata); - KapaPlotVector (kapa, nplot, dXplot, "x"); - KapaPlotVector (kapa, nplot, dYplot, "y"); - - //Overplot bounding box, peak of distribution - float xbound[5] = { -maxOffpix, maxOffpix, maxOffpix, -maxOffpix, -maxOffpix}; - float ybound[5] = { -maxOffpix, -maxOffpix, maxOffpix, maxOffpix, -maxOffpix}; - float xbin[5] = {stats->offset.x - 0.5 * Scale, stats->offset.x + 0.5 * Scale, stats->offset.x + 0.5 * Scale, - stats->offset.x - 0.5 * Scale, stats->offset.x - 0.5 * Scale}; - float ybin[5] = {stats->offset.y - 0.5 * Scale, stats->offset.y - 0.5 * Scale, stats->offset.y + 0.5 * Scale, - stats->offset.y + 0.5 * Scale, stats->offset.y - 0.5 * Scale}; - graphdata.color = KapaColorByName("red"); - graphdata.style = 0; - graphdata.size = 1.0; - KapaPrepPlot(kapa, 5, &graphdata); - KapaPlotVector (kapa, 5, xbound, "x"); - KapaPlotVector (kapa, 5, ybound, "y"); - KapaPrepPlot(kapa, 5, &graphdata); - KapaPlotVector (kapa, 5, xbin, "x"); - KapaPlotVector (kapa, 5, ybin, "y"); - - KapaPNG(kapa, "dXdY.png"); - char c; - fscanf(stdin, "%c", &c); - /*** end ***/ -#endif psFree (sort); @@ -767,5 +704,4 @@ PS_ASSERT_PTR_NON_NULL(ref, NULL); PS_ASSERT_PTR_NON_NULL(config, NULL); - bool status; double xMin, xMax, yMin, yMax; Index: /branches/cnb_branch_20081011/psModules/src/astrom/pmAstrometryObjects.h =================================================================== --- /branches/cnb_branch_20081011/psModules/src/astrom/pmAstrometryObjects.h (revision 20334) +++ /branches/cnb_branch_20081011/psModules/src/astrom/pmAstrometryObjects.h (revision 20335) @@ -4,6 +4,6 @@ * @author EAM, IfA * - * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-11-21 07:02:55 $ + * @version $Revision: 1.17.42.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-10-22 23:35:13 $ * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii */ @@ -238,11 +238,9 @@ * */ -bool pmAstromFitFPA( - pmFPA *fpa, - psArray *st1, - psArray *st2, - psArray *match, - psMetadata *config -); +bool pmAstromFitFPA(pmFPA *fpa, + psArray *st1, + psArray *st2, + psArray *match, + psMetadata *config); @@ -269,6 +267,5 @@ psArray *st2, psArray *match, - psMetadata *config -); + psMetadata *config); Index: /branches/cnb_branch_20081011/psModules/src/psmodules.h =================================================================== --- /branches/cnb_branch_20081011/psModules/src/psmodules.h (revision 20334) +++ /branches/cnb_branch_20081011/psModules/src/psmodules.h (revision 20335) @@ -81,4 +81,5 @@ #include #include +#include // the following headers are from psModule:imcombine Index: /branches/cnb_branch_20081011/psastro/src/psastro.h =================================================================== --- /branches/cnb_branch_20081011/psastro/src/psastro.h (revision 20334) +++ /branches/cnb_branch_20081011/psastro/src/psastro.h (revision 20335) @@ -51,4 +51,5 @@ psArray *psastroRemoveClumps (psArray *input, int scale); + // utility functions: bool psastroUpdateChipToFPA (pmFPA *fpa, pmChip *chip, psArray *rawstars, psArray *refstars); @@ -84,4 +85,6 @@ 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 psastroVisualPlotMosaicMatches (psArray *rawstars, psArray *refstars, psArray *match, int iteration, psMetadata *recipe); // demo plots Index: /branches/cnb_branch_20081011/psastro/src/psastroMosaicOneChip.c =================================================================== --- /branches/cnb_branch_20081011/psastro/src/psastroMosaicOneChip.c (revision 20334) +++ /branches/cnb_branch_20081011/psastro/src/psastroMosaicOneChip.c (revision 20335) @@ -5,5 +5,5 @@ if (!status) { \ psError(PSASTRO_ERR_CONFIG, false, MESSAGE); \ - return false; } + return false; } bool psastroMosaicOneChip (pmChip *chip, pmReadout *readout, psMetadata *recipe, psMetadata *updates, int iteration) { @@ -29,5 +29,5 @@ // correct radius to FP units (physical pixel scale in microns per pixel) - REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32, "Failed to lookup pixel scale"); + REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32, "Failed to lookup pixel scale"); // allowed limits for valid solutions @@ -42,5 +42,5 @@ int order = psMetadataLookupS32 (&status, recipe, orderWord); if (!status || (order == -1)) { - order = defaultOrder; + order = defaultOrder; } @@ -50,7 +50,7 @@ if ((match->n < 8) && (order >= 1)) order = 0; if ((match->n < 3) || (order < 0) || (order > 3)) { - psLogMsg ("psastro", 3, "insufficient stars (%ld) or invalid order (%d)", match->n, order); - return false; - } + psLogMsg ("psastro", 3, "insufficient stars (%ld) or invalid order (%d)", match->n, order); + return false; + } psLogMsg ("psastro", PS_LOG_DETAIL, "mosaic fit chip order %d", order); @@ -60,25 +60,25 @@ // coefficients frozen to the current values if (order == 0) { - // set FIT mask for all higher order terms of the existing solution - // any existing SET masks will be retained. - for (int i = 0; i <= chip->toFPA->x->nX; i++) { - for (int j = 0; j <= chip->toFPA->x->nY; j++) { - if (i + j > 0) { - chip->toFPA->x->coeffMask[i][j] |= PS_POLY_MASK_FIT; - chip->toFPA->y->coeffMask[i][j] |= PS_POLY_MASK_FIT; - } - } - } + // set FIT mask for all higher order terms of the existing solution + // any existing SET masks will be retained. + for (int i = 0; i <= chip->toFPA->x->nX; i++) { + for (int j = 0; j <= chip->toFPA->x->nY; j++) { + if (i + j > 0) { + chip->toFPA->x->coeffMask[i][j] |= PS_POLY_MASK_FIT; + chip->toFPA->y->coeffMask[i][j] |= PS_POLY_MASK_FIT; + } + } + } } else { - 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; - } - } - } + 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; + } + } + } } @@ -87,11 +87,11 @@ psStats *fitStats = NULL; // if (order == 0) { -// fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); -// fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA"); -// fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER"); +// fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); +// fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA"); +// fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER"); // } else { - fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); - fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA"); - fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER"); + fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); + fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA"); + fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER"); // } @@ -99,6 +99,6 @@ pmAstromFitResults *results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats); if (!results) { - psError(PSASTRO_ERR_DATA, false, "failed to perform the matched fit\n"); - return false; + psError(PSASTRO_ERR_DATA, false, "failed to perform the matched fit\n"); + return false; } @@ -124,10 +124,10 @@ psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d", astError, astNstar); if ((maxError > 0) && (astError > maxError)) { - psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError); - validSolution = false; + psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError); + validSolution = false; } if (astNstar < minNstar) { - psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar); - validSolution = false; + psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar); + validSolution = false; } @@ -136,9 +136,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, "", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars @@ -147,4 +147,7 @@ psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); + //plot results + psastroVisualPlotMosaicOneChip(rawstars, refstars, match, recipe); + psFree (fitStats); psFree (results); Index: /branches/cnb_branch_20081011/psastro/src/psastroMosaicSetMatch.c =================================================================== --- /branches/cnb_branch_20081011/psastro/src/psastroMosaicSetMatch.c (revision 20334) +++ /branches/cnb_branch_20081011/psastro/src/psastroMosaicSetMatch.c (revision 20335) @@ -10,22 +10,22 @@ // use small radius to match stars (assume starting astrometry is good) - bool status = false; + bool status = false; sprintf (radiusWord, "PSASTRO.MOSAIC.RADIUS.N%d", iteration); - double RADIUS = psMetadataLookupF32 (&status, recipe, radiusWord); - if (!status) { - psError(PS_ERR_IO, false, "Failed to lookup matching radius: %s", radiusWord); - psFree (view); - return false; - } + double RADIUS = psMetadataLookupF32 (&status, recipe, radiusWord); + if (!status) { + psError(PS_ERR_IO, false, "Failed to lookup matching radius: %s", radiusWord); + psFree (view); + return false; + } if (RADIUS <= 0.0) { - if (iteration == 0) { - psError(PS_ERR_IO, false, "Invalid match radius for first iteration: %s", radiusWord); - psFree (view); - return false; - } - psWarning ("skipping match for iteration %d\n", iteration); - psFree (view); - return true; + if (iteration == 0) { + psError(PS_ERR_IO, false, "Invalid match radius for first iteration: %s", radiusWord); + psFree (view); + return false; + } + psWarning ("skipping match for iteration %d\n", iteration); + psFree (view); + return true; } @@ -34,32 +34,34 @@ psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); if (!chip->process || !chip->file_exists) { continue; } - if (!chip->fromFPA) { continue; } - - while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { + if (!chip->fromFPA) { continue; } + + while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); if (!cell->process || !cell->file_exists) { continue; } - // process each of the readouts - // XXX there can only be one readout per chip, right? - while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { - if (! readout->data_exists) { continue; } + // process each of the readouts + // XXX there can only be one readout per chip, right? + while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { + if (! readout->data_exists) { continue; } - // 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 *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; } - psTrace ("psastro", 4, "Trying %ld refstars\n", refstars->n); + // select the raw objects for this readout + psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); + if (refstars == NULL) { continue; } + psTrace ("psastro", 4, "Trying %ld refstars\n", refstars->n); - psArray *matches = pmAstromRadiusMatchChip (rawstars, refstars, RADIUS); - psTrace ("psastro", 4, "Matched %ld refstars\n", matches->n); + psArray *matches = pmAstromRadiusMatchChip (rawstars, refstars, RADIUS); + psTrace ("psastro", 4, "Matched %ld refstars\n", matches->n); - // XXX drop the old one - psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.MATCH", PS_DATA_ARRAY | PS_META_REPLACE, "astrometry matches", matches); - psFree (matches); - } - } + psastroVisualPlotMosaicMatches(rawstars, refstars, matches, iteration, recipe); + + // XXX drop the old one + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.MATCH", PS_DATA_ARRAY | PS_META_REPLACE, "astrometry matches", matches); + psFree (matches); + } + } } psFree (view); Index: /branches/cnb_branch_20081011/psastro/src/psastroVisual.c =================================================================== --- /branches/cnb_branch_20081011/psastro/src/psastroVisual.c (revision 20334) +++ /branches/cnb_branch_20081011/psastro/src/psastroVisual.c (revision 20335) @@ -1,5 +1,9 @@ -/***Diagnostic plots for psastro***/ - +/***********************************/ +/***Diagnostic plots for psastro****/ +/***********************************/ # include "psastroInternal.h" +# include +# include +# include # if (HAVE_KAPA) @@ -11,14 +15,16 @@ //variables to determine when things are plotted static bool isVisual = false; -static bool plotRawStars = false; -static bool plotRefStars = false; -static bool plotLumFunc = false; +static bool plotRawStars = true; +static bool plotRefStars = true; +static bool plotLumFunc = true; static bool plotRemoveClumps = true; -static bool plotOneChipFit = false; -static bool plotAstromGuessCheck = false; +static bool plotOneChipFit = true; +static bool plotAstromGuessCheck = true; +static bool plotMosaicMatches = true; +static bool plotMosaicOneChip = true; // variables to store plotting window indices -static int kapa = -1; //for plots -static int kapa2 = -1; //also for plots +static int kapa = -1; +static int kapa2 = -1; // helper routine prototypes @@ -27,12 +33,12 @@ bool psastroVisualTripleOverplot (int kapa, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing); bool psastroVisualCreateScaleVec (psVector *zVec, psVector *zScale, bool increasing); - - +bool psastroVisualResidPlot (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe, char *title); + + +/*****************************/ /** Initialization Routines **/ - - -/** - * start or stop plotting - */ +/*****************************/ + +/** start or stop plotting */ bool psastroSetVisual (bool mode) { isVisual = mode; @@ -41,7 +47,5 @@ -/** - * open, name, and resize a window if necessary - */ +/** open, name, and resize a window if necessary*/ bool psastroVisualInitWindow( int *kapid, char *name ) { if (*kapid == -1) { @@ -58,7 +62,5 @@ -/** - * ask the user how to proceed - */ +/** ask the user how to proceed*/ bool psastroVisualAskUser( bool *plotflag ) { char key[10]; @@ -77,7 +79,5 @@ -/** - * destroy windows at the end of a run - */ +/** destroy windows at the end of a run*/ bool psastroVisualClose() { if(kapa != -1) @@ -88,22 +88,24 @@ } + +/***********************/ /** Plotting Routines **/ - -/* psastroVisualPlotRawStars */ -// Show raw star astrometry as determined from first pass fit -// Invoked from within psastroAstromGuess +/***********************/ + + +/** psastroVisualPlotRawStars + * Plot raw stars as determined from first pass astrometry fit + * Called within psastroAstromGeuss + */ bool psastroVisualPlotRawStars (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 (!psastroVisualInitWindow (&kapa, "psastro:plots")) return false; - Graphdata graphdata; KapaSection section; - - bool status = false; - float iMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MIN"); - float iMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MAX"); KapaInitGraph (&graphdata); @@ -118,4 +120,9 @@ 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); @@ -128,4 +135,6 @@ KapaSetSection (kapa, §ion); psFree (section.name); + + //Chip coordinates int n = 0; for (int i = 0; i < rawstars->n; i++) { @@ -147,4 +156,5 @@ psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false); + //Focal Plane Coordinates section.x = 0.5; section.y = 0.5; @@ -173,4 +183,5 @@ psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false); + // Tangent Plane Coordinates section.x = 0.0; section.y = 0.0; @@ -199,4 +210,5 @@ psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false); + //sky coordinates section.x = 0.5; section.y = 0.0; @@ -250,17 +262,17 @@ -/* psastroVisualPlotRefStars */ -// plot the location of references stars over the entire fpa -// invoked during psastroChooseRefStars +/** psastroVisualPlotRefStars + * plot the location of references stars over the entire fpa + * invoked during psastroChooseRefStars + */ bool psastroVisualPlotRefStars (psArray *refstars, psMetadata *recipe) { + //make sure we want to plot this if (!isVisual || !plotRefStars) return true; Graphdata graphdata; + + //set up plotting variables if (!psastroVisualInitWindow (&kapa, "psastro:plots")) return false; - - bool status = false; - float rMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MIN"); - float rMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MAX"); KapaInitGraph (&graphdata); @@ -271,4 +283,9 @@ 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); @@ -309,17 +326,19 @@ -/* psastroVisualPlotLuminosityFunction */ -// Plot the two luminosity functions created. -// Called from psastroRefStarSubset. -// The luminosity functions are used to select a subset of reference stars, -// so we plot the cutoff that defines this subset +/** psastroVisualPlotLuminosityFunction + * 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 psastroVisualPlotLuminosityFunction (psVector *lnMag, psVector *Mag, pmLumFunc *lumFunc, pmLumFunc *rawFunc) { + // make sure we want to plot this if ( !isVisual || !plotLumFunc ) return true; + + //set up plotting variables if ( !psastroVisualInitWindow (&kapa, "psastro:plots")) return false; - //Set up graph window Graphdata graphdata; KapaSection section1 = {"s1", .1, .1, .85, .35}; @@ -396,18 +415,19 @@ -/** +/** psastroVisualPlotRemoveClumps * 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 psastroVisualPlotRemoveClumps (psArray *input, psImage *count, int scale, float limit) { + + //make sure we want to plot this if (!isVisual || !plotRemoveClumps) return true; + + //set up plot variables if ( !psastroVisualInitWindow (&kapa, "psastro:plots")) return false; KapaSection section; Graphdata graphdata; - - //initialize graph information section.x = 0; section.dx = .9; @@ -419,9 +439,10 @@ 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 @@ -486,4 +507,6 @@ } } + + //ask for user input and finish psastroVisualAskUser (&plotRemoveClumps); psFree (xVec); @@ -494,15 +517,353 @@ -/* psastroVisualPlotOneChipFit */ -// assess the goodness of fit for a signle chip by -// plotting the fit residuals -// invoked during psastroOneChipFit +/** psastroVisualPlotOneChipFit + * assess the goodness of fit for a signle chip by + * plotting the fit residuals + * invoked during psastroOneChipFit + */ bool psastroVisualPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe) { + + //make sure we want to plot this if (!isVisual || !plotOneChipFit) return true; - if (!psastroVisualInitWindow (&kapa, "psastro:plots")) + //plot the residuals + if (!psastroVisualResidPlot(rawstars, refstars, match, recipe, "Single Chip Fit Residuals (Chip Coordinates)")) return false; + //ask for user input and finish + psastroVisualAskUser(&plotOneChipFit); + return true; +} + + +/** + * 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 psastroVisualPlotAstromGuessCheck (psVector *cornerPo, psVector *cornerQo, + psVector *cornerPn, psVector *cornerQn, + psVector *cornerPd, psVector *cornerQd) { + + //make sure we want to plot this + if (!isVisual || !plotAstromGuessCheck) return true; + + //set up graph window + if ( !psastroVisualInitWindow (&kapa, "psastro:plots")) + 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); + + psastroVisualScaleGraphdata (&graphdata, cornerPo, cornerPo); + 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); + + psastroVisualScaleGraphdata (&graphdata, xResid, yResid); + 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); + psastroVisualAskUser (&plotAstromGuessCheck); + return true; +} + + +/** + * psastroVisualPlotMosaicOneChip (rawstars, refstars, match); + * plot the residuals between raw stars and ref stars after + * fitting in psastroMosaicOneChip + */ +bool psastroVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars, + psArray *match, psMetadata *recipe) +{ + + //make sure we want to plot this + if (!isVisual || !plotMosaicOneChip) return false; + + //plot the residuals + if (!psastroVisualResidPlot(rawstars, refstars, match, recipe, "Single Chip Fit Residuals - Mosaic Mode")) + return false; + + //ask for user input and finish + psastroVisualAskUser(&plotMosaicOneChip); + + return true; +} + +/** psastroVisualPlotMosaicMatches + * Plot the matches between raw and reference stars during psastroVisualMosaicSetMatch + */ +bool psastroVisualPlotMosaicMatches( psArray *rawstars, psArray *refstars, + psArray *match, int iteration, + psMetadata *recipe) +{ + //make sure we want to plot this + if (!isVisual || !plotMosaicMatches) return true; + char *title; + switch(iteration) { + case 0: + title = "Matches found during psastroMosaicSetMatch iteration 0"; + break; + case 1 : + title = "Matches found during psastroMosaicSetMatch iteration 1"; + break; + case 2: + title = "Matches found during psastroMosaicSetMatch iteration 2"; + break; + case 3: + title = "Matches found during psastroMosaicSetMatch iteration 3"; + break; + case 4: + title = "Matches found during psastroMosaicSetMatch iteration 4"; + break; + default: + title = "Matches found during psastroMosaicSetMatch"; + break; + } + + if (!psastroVisualResidPlot(rawstars, refstars, match, recipe, title)) + return false; + + //ask for user input + psastroVisualAskUser (&plotMosaicMatches); + return true; +} + + +/*********************/ +/** Helper Routines **/ +/*********************/ + + +/** psastroVisualTriplePlot + * plot 2 vectors whose point sizes are scaled by a third vector + * autoscales the plot + */ +bool psastroVisualTriplePlot (int kapid, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing) +{ + psastroVisualScaleGraphdata (graphdata, xVec, yVec); + //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); + psastroVisualCreateScaleVec (zVec, zScale, increasing); + + KapaSetFont (kapid, "helvetica", 14); + KapaSetLimits(kapid, graphdata); + KapaBox (kapid, graphdata); + + // 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; +} + + +/** psastroVisualTripleOverplot + * same as above, but does not rescale the plot or redraw the bounding box + */ +bool psastroVisualTripleOverplot (int kapid, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing) { + + // set the scale vector + psVector *zScale = psVectorAlloc (zVec->n, PS_DATA_F32); + psastroVisualCreateScaleVec (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; +} + + +/** psastroVisualCreateScaleVec + * create a vector of plot point sizes from a vector + */ +bool psastroVisualCreateScaleVec (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; +} + + +/** psastroVisualScaleGraphdata + * Scale the graphdata structure based on x and y coordinates. Use sigma clipping to + * prevent outliers from making te plot region too big. + */ +bool psastroVisualScaleGraphdata(Graphdata *graphdata, psVector *xVec, psVector *yVec) { + + 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; + + // 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); + return true; +} + + +/** psastroVisualResdiPlot + * Plots the residuals between matched raw and reference stars. Used to assess the quality of an + * astrometry solution + */ +bool psastroVisualResidPlot (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe, char *title) { + + //set up the first window + if (!psastroVisualInitWindow( &kapa, "psastro:plots")) return false; + + //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"); @@ -511,19 +872,4 @@ float rMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MAX"); - //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; - psVector *xVec = psVectorAlloc (match->n, PS_TYPE_F32); psVector *yVec = psVectorAlloc (match->n, PS_TYPE_F32); @@ -661,9 +1007,9 @@ psStringAppend (§ion.name, "a5"); KapaSetSection (kapa, §ion); - KapaSendLabel (kapa, "Single Chip Fit Residuals (Chip Coordinates)", KAPA_LABEL_XP); + KapaSendLabel (kapa, title, KAPA_LABEL_XP); psFree (section.name); - // *** X vs Y plot (different window) + // X vs Y plot (different window) if (!psastroVisualInitWindow( &kapa2, "psastro:plots")) return false; @@ -774,6 +1120,4 @@ psastroVisualTripleOverplot (kapa2, &graphdata, xVec, yVec, zVec, false); - psastroVisualAskUser (&plotOneChipFit); - psFree (xVec); psFree (yVec); @@ -782,237 +1126,4 @@ } - -/** - * 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 psastroVisualPlotAstromGuessCheck (psVector *cornerPo, psVector *cornerQo, - psVector *cornerPn, psVector *cornerQn, - psVector *cornerPd, psVector *cornerQd) { - - if (!isVisual || !plotAstromGuessCheck) return true; - if ( !psastroVisualInitWindow (&kapa, "psastro:plots")) - return false; - - //set up graph window - 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); - - psastroVisualScaleGraphdata (&graphdata, cornerPo, cornerPo); - 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] / cornerPn->data.F32[i]); - yResid->data.F32[i] = (cornerQd->data.F32[i] / cornerQn->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); - - psastroVisualScaleGraphdata (&graphdata, xResid, yResid); - KapaSetLimits (kapa, &graphdata); - KapaBox (kapa, &graphdata); - KapaSendLabel (kapa, "dP / P", KAPA_LABEL_XM); - KapaSendLabel (kapa, "dQ / Q", 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); - psastroVisualAskUser (&plotAstromGuessCheck); - return true; -} - -/** Helper Routines **/ - - -/* psastroVisualTriplePlot */ -// plot 2 vectors whose point sizes are scaled by a third vector -// autoscales the plot -bool psastroVisualTriplePlot (int kapid, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing) -{ - psastroVisualScaleGraphdata (graphdata, xVec, yVec); - - // set the scale vector - psVector *zScale = psVectorAlloc (zVec->n, PS_DATA_F32); - psastroVisualCreateScaleVec (zVec, zScale, increasing); - - KapaSetFont (kapid, "helvetica", 14); - KapaSetLimits(kapid, graphdata); - KapaBox (kapid, graphdata); - - // 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; -} - - -/* psastroVisualTripleOverplot */ -// same as above, but does not rescale the plot or redraw the bounding box -bool psastroVisualTripleOverplot (int kapid, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing) { - - // set the scale vector - psVector *zScale = psVectorAlloc (zVec->n, PS_DATA_F32); - psastroVisualCreateScaleVec (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; -} - - -/* psastroVisualCreateScaleVec */ -//create a vector of plot point sizes from a vector -bool psastroVisualCreateScaleVec (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; -} - - -/** psastroVisualScaleGraphdata - * Scale the graphdata structure based on x and y coordinates. Use sigma clipping to - * prevent outliers from making te plot region too big. - */ -bool psastroVisualScaleGraphdata(Graphdata *graphdata, psVector *xVec, psVector *yVec) { - graphdata->xmin = +FLT_MAX; - graphdata->xmax = -FLT_MAX; - graphdata->ymin = +FLT_MAX; - graphdata->ymax = -FLT_MAX; - bool baddataX = true; - bool baddataY = true; - - //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; - - 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; - baddataX = false; - 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(yVec->data.F32[i])) continue; - if (yVec->data.F32[i] > yhi || yVec->data.F32[i] < ylo) continue; - baddataY = false; - graphdata->ymin = PS_MIN (graphdata->ymin, yVec->data.F32[i]); - graphdata->ymax = PS_MAX (graphdata->ymax, yVec->data.F32[i]); - } - - // check to see if there was no good data - if ( baddataX || baddataY) { - graphdata->xmin = -1; - graphdata->ymin = -1; - graphdata->xmax = 1; - graphdata->ymax = 1; - return false; - } - 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; -} - # endif