Index: trunk/psModules/src/extras/.cvsignore =================================================================== --- trunk/psModules/src/extras/.cvsignore (revision 13896) +++ 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 13896) +++ 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