Index: /trunk/psphot/src/Makefile.am =================================================================== --- /trunk/psphot/src/Makefile.am (revision 10392) +++ /trunk/psphot/src/Makefile.am (revision 10393) @@ -3,19 +3,10 @@ libpsphot_la_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSPHOT_CFLAGS) -# bin_PROGRAMS = psphot polyfitTest growthTest bin_PROGRAMS = psphot psphot_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSPHOT_CFLAGS) -psphot_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) +psphot_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PSPHOT_LIBS) psphot_LDADD = libpsphot.la psphot_SOURCES = psphot.c - -polyfitTest_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSPHOT_CFLAGS) -polyfitTest_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) -polyfitTest_SOURCES = polyfitTest.c - -growthTest_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSPHOT_CFLAGS) -growthTest_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) -growthTest_SOURCES = growthTest.c psphotGrowthCurve.c libpsphot_la_SOURCES = \ @@ -53,5 +44,6 @@ psphotWeightBias.c \ psphotSourceFreePixels.c \ - psphotTestPSF.c \ + psphotSummaryPlots.c \ + psphotTestPSF.c \ psphotCleanup.c Index: unk/psphot/src/growthTest.c =================================================================== --- /trunk/psphot/src/growthTest.c (revision 10392) +++ (revision ) @@ -1,89 +1,0 @@ -# include "psphot.h" - -int main (int argc, char **argv) { - - if (argc != 3) { - fprintf (stderr, "USAGE: growthTest\n"); - exit (2); - } - - // does the growth curve correctly fix aperture mags? - pmModelGroupInit (); - - // generate a simple readout - pmReadout *readout = pmReadoutAlloc (NULL); - readout->image = psImageAlloc (512, 512, PS_TYPE_F32); - readout->mask = psImageAlloc (512, 512, PS_TYPE_U8); - - // generate a simple psf - pmModelType type = pmModelSetType ("PS_MODEL_GAUSS"); - psPolynomial2D *psfTrend = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 0, 0); - pmPSF *psf = pmPSFAlloc (type, true, psfTrend); - psf->growth = pmGrowthCurveAlloc (2.0, 100.0, 15.0); - - psPolynomial2D *poly; - - poly = psf->params_NEW->data[PM_PAR_SXX]; - poly->coeff[0][0] = M_SQRT2*1.5; - - poly = psf->params_NEW->data[PM_PAR_SYY]; - poly->coeff[0][0] = M_SQRT2*1.5; - - poly = psf->params_NEW->data[PM_PAR_SXY]; - poly->coeff[0][0] = 0; - - // create template model - pmModel *modelRef = pmModelAlloc(psf->type); - modelRef->params->data.F32[PM_PAR_SKY] = 0; - modelRef->params->data.F32[PM_PAR_I0] = 1000; - - // generate the growth curve (element of psf) - psphotGrowthCurve (readout, psf, false); - - for (float radius = 2; radius < 8; radius += 1.0) { - for (float dx = -1; dx <= 1; dx += 0.1) { - for (float dy = -1; dy <= 1; dy += 0.1) { - - float fitMag, apMag, deltaMag; - - // assign the x and y coords to the image center - // create an object with center intensity of 1000 - modelRef->params->data.F32[PM_PAR_XPOS] = 200+dx; - modelRef->params->data.F32[PM_PAR_YPOS] = 200+dy; - - // create modelPSF from this model - pmModel *model = pmModelFromPSF (modelRef, psf); - model->radiusFit = radius; - - // measure the fitMag for this model - pmSourcePhotometryModel (&fitMag, model); - - psRegion region = {180, 220, 180, 220}; - psImage *imageView = psImageSubset (readout->image, region); - psImage *maskView = psImageSubset (readout->mask, region); - - psImageInit (imageView, 0.0); - psImageInit (maskView, 0); - pmModelAdd (imageView, maskView, model, false, false); - - float dX = 0.5 - model->params->data.F32[PM_PAR_XPOS] + (int)model->params->data.F32[PM_PAR_XPOS]; - float dY = 0.5 - model->params->data.F32[PM_PAR_YPOS] + (int)model->params->data.F32[PM_PAR_YPOS]; - psImage *tmpImage = psImageShift (NULL, imageView, dX, dY, 0.0, PS_INTERPOLATE_BICUBE); - // psImage *tmpImage = tmp_psImageShift (imageView, dX, dY); - - psImageKeepCircle (maskView, 200+dx+dX, 200+dy+dY, radius, "OR", PM_MASK_MARK); - pmSourcePhotometryAper (&apMag, model, tmpImage, maskView); - psImageKeepCircle (maskView, 200+dx+dX, 200+dy+dY, radius, "AND", PS_NOT_U8(PM_MASK_MARK)); - - deltaMag = pmGrowthCurveCorrect (psf->growth, radius); - - fprintf (stdout, "%f %f %f %f %f %f\n", radius, dx, dy, fitMag, apMag, deltaMag); - - psFree (tmpImage); - psFree (imageView); - psFree (maskView); - psFree (model); - } - } - } -} Index: unk/psphot/src/polyfitTest.c =================================================================== --- /trunk/psphot/src/polyfitTest.c (revision 10392) +++ (revision ) @@ -1,93 +1,0 @@ -# include "psphot.h" - -int main (int argc, char **argv) { - - float X, Y, SX, SY, SXY; - - if (argc != 3) { - fprintf (stderr, "USAGE: polyfitTest (input) (output)\n"); - exit (2); - } - - FILE *f = fopen (argv[1], "r"); - if (f == NULL) exit (1); - - psVector *x = psVectorAlloc (100, PS_TYPE_F32); - psVector *y = psVectorAlloc (100, PS_TYPE_F32); - psVector *sx = psVectorAlloc (100, PS_TYPE_F32); - psVector *sy = psVectorAlloc (100, PS_TYPE_F32); - psVector *sxy = psVectorAlloc (100, PS_TYPE_F32); - - int i = 0; - while (fscanf (f, "%f %f %f %f %f %*f %*d", &X, &Y, &SX, &SY, &SXY) != EOF) { - x->data.F32[i] = X; - y->data.F32[i] = Y; - sx->data.F32[i] = SX; - sy->data.F32[i] = SY; - sxy->data.F32[i] = SXY; - - psVectorExtend (x, 100, 1); - psVectorExtend (y, 100, 1); - psVectorExtend (sx, 100, 1); - psVectorExtend (sy, 100, 1); - psVectorExtend (sxy, 100, 1); - - i++; - } - fclose (f); - - fprintf (stderr, "loaded %ld pts\n", x->n); - - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); - psVector *mask = psVectorAlloc (x->n, PS_TYPE_MASK); - mask->n = x->n; - psVectorInit (mask, 0); - - psPolynomial2D *polySX = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2); - polySX->mask[2][1] = 1; - polySX->mask[2][2] = 1; - polySX->mask[1][2] = 1; - - psVectorClipFitPolynomial2D (polySX, stats, mask, 0xff, sx, NULL, x, y); - fprintf (stderr, "stats: %f +/- %f\n", stats->robustMedian, stats->robustStdev); - fprintf (stderr, "SX: %g %g %g\n", polySX->coeff[0][0], polySX->coeff[1][0], polySX->coeff[0][1]); - fprintf (stderr, "SX: %g %g %g\n", polySX->coeff[0][2], polySX->coeff[1][1], polySX->coeff[2][0]); - - psPolynomial2D *polySY = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2); - polySY->mask[2][1] = 1; - polySY->mask[2][2] = 1; - polySY->mask[1][2] = 1; - - psVectorClipFitPolynomial2D (polySY, stats, mask, 0xff, sy, NULL, x, y); - fprintf (stderr, "stats: %f +/- %f\n", stats->robustMedian, stats->robustStdev); - fprintf (stderr, "SY: %g %g %g\n", polySY->coeff[0][0], polySY->coeff[1][0], polySY->coeff[0][1]); - fprintf (stderr, "SY: %g %g %g\n", polySY->coeff[0][2], polySY->coeff[1][1], polySY->coeff[2][0]); - - psPolynomial2D *polySXY = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2); - polySXY->mask[2][1] = 1; - polySXY->mask[2][2] = 1; - polySXY->mask[1][2] = 1; - - psVectorClipFitPolynomial2D (polySXY, stats, mask, 0xff, sxy, NULL, x, y); - fprintf (stderr, "stats: %f +/- %f\n", stats->robustMedian, stats->robustStdev); - fprintf (stderr, "SXY: %g %g %g\n", polySXY->coeff[0][0], polySXY->coeff[1][0], polySXY->coeff[0][1]); - fprintf (stderr, "SXY: %g %g %g\n", polySXY->coeff[0][2], polySXY->coeff[1][1], polySXY->coeff[2][0]); - - psVector *sxf = psPolynomial2DEvalVector (polySX, x, y); - psVector *syf = psPolynomial2DEvalVector (polySY, x, y); - psVector *sxyf = psPolynomial2DEvalVector (polySXY, x, y); - - f = fopen (argv[2], "w"); - for (i = 0; i < x->n; i++) { - fprintf (f, "%f %f %f %f %f %f %f %f\n", - x->data.F32[i], y->data.F32[i], - sx->data.F32[i], sy->data.F32[i], sxy->data.F32[i], - sxf->data.F32[i], syf->data.F32[i], sxyf->data.F32[i]); - } - fclose (f); - exit (0); -} - - // 0,0 1,0 2,0 - // 0,1 1,1 2,1 - // 0,2 1,2 2,2 Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 10392) +++ /trunk/psphot/src/psphot.h (revision 10393) @@ -92,2 +92,9 @@ bool psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, pmSourceFitMode mode); +// plotting functions (available if libkapa is installed) +bool psphotPlotMoments (pmConfig *config, pmFPAview *view, psArray *sources); +bool psphotPlotPSFModel (pmConfig *config, pmFPAview *view, psArray *sources); + +// move to psLib or psModules +int psphotKapaOpen (); +bool psphotKapaClose (); Index: /trunk/psphot/src/psphotParseCamera.c =================================================================== --- /trunk/psphot/src/psphotParseCamera.c (revision 10392) +++ /trunk/psphot/src/psphotParseCamera.c (revision 10393) @@ -66,4 +66,10 @@ } + // optionally save the PSF Model + if (psMetadataLookupBool(NULL, recipe, "SAVE.PLOTS")) { + pmFPAfileDefineOutput (config, input->fpa, "PSPHOT.MOMENT.PLT"); + pmFPAfileDefineOutput (config, input->fpa, "PSPHOT.PSFMODEL.PLT"); + } + // pmFPAfileConstruct (config->files, format, config->camera, "PSPHOT.PSF_SAMPLE"); Index: /trunk/psphot/src/psphotReadout.c =================================================================== --- /trunk/psphot/src/psphotReadout.c (revision 10392) +++ /trunk/psphot/src/psphotReadout.c (revision 10393) @@ -70,4 +70,7 @@ return false; } + // XXX test the plotting function + psphotPlotMoments (config, view, sources); + if (!strcasecmp (breakPt, "MOMENTS")) { if (!psphotReadoutCleanup (config, readout, recipe, NULL, sources)) return false; @@ -117,4 +120,7 @@ // linear PSF fit to remaining peaks psphotEnsemblePSF (readout, sources, recipe, psf, TRUE); + + // XXX test the plotting function + psphotPlotPSFModel (config, view, sources); finish: @@ -169,4 +175,6 @@ pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND"); + psphotKapaClose (); + psFree (psf); psFree (header); Index: /trunk/psphot/src/psphotSummaryPlots.c =================================================================== --- /trunk/psphot/src/psphotSummaryPlots.c (revision 10393) +++ /trunk/psphot/src/psphotSummaryPlots.c (revision 10393) @@ -0,0 +1,299 @@ +# include "psphot.h" + +// the top portion of this file defines plotting functions which use kapa for plotting. +// if kapa is not available, these functions are defined in the bottom portion as stubs +// which perform NOP and return false (XXX should this be true??) + +# if (HAVE_KAPA) + +# include + +// XXX not thread safe (perhaps not needed) +// to make this thread safe, we could check the thread ID and tie to it +static int kapa_fd = -1; +int psphotKapaOpen () { + + if (kapa_fd == -1) { + // XXX make -noX an option + // XXX -noX is crashing kapa on point plotting + kapa_fd = KapaOpen ("kapa", "psphot"); + } + return kapa_fd; +} + +bool psphotKapaClose () { + + if (kapa_fd == -1) return true; + KapaClose (kapa_fd); + kapa_fd = -1; + return true; +} + +// plot the sx, sy moments plane (faint and bright sources) +bool psphotPlotMoments (pmConfig *config, pmFPAview *view, psArray *sources) { + + // select model pixels (from output background model file, or create internal file) + pmFPAfile *file = psMetadataLookupPtr (NULL, config->files, "PSPHOT.MOMENT.PLT"); + if (file == NULL) { + psLogMsg ("psphot", 3, "skipping moments plot"); + return false; + } + + // pmFPAfileOpen defers disk I/O for KAPA files: just get the correct name + pmFPAfileOpen (file, view, config); + + Graphdata graphdata; + + psLogMsg ("psphot", 3, "creating moments plot"); + + int kapa = psphotKapaOpen (); + if (kapa == -1) { + psError(PSPHOT_ERR_UNKNOWN, true, "failure to open kapa"); + return false; + } + + KapaResize (kapa, 500, 500); + KapaInitGraph (&graphdata); + + // examine sources to set data range + graphdata.xmin = -0.05; + graphdata.ymin = -0.05; + graphdata.xmax = +2.05; + graphdata.ymax = +2.05; + KapaSetLimits (kapa, &graphdata); + + KapaSetFont (kapa, "helvetica", 14); + KapaBox (kapa, &graphdata); + KapaSendLabel (kapa, "&ss&h_x| (pixels)", KAPA_LABEL_XM); + KapaSendLabel (kapa, "&ss&h_y| (pixels)", KAPA_LABEL_YM); + + psVector *xBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *yBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *xFaint = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *yFaint = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + + // construct the vectors + int nB = 0; + int nF = 0; + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (source->moments == NULL) continue; + + xFaint->data.F32[nF] = source->moments->Sx; + yFaint->data.F32[nF] = source->moments->Sy; + nF++; + + // XXX make this a user-defined cutoff + if (source->moments->SN < 25) continue; + + xBright->data.F32[nB] = source->moments->Sx; + yBright->data.F32[nB] = source->moments->Sy; + nB++; + } + xFaint->n = nF; + yFaint->n = nF; + + xBright->n = nB; + yBright->n = nB; + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 0; + graphdata.size = 0.3; + graphdata.style = 2; + KapaPrepPlot (kapa, nF, &graphdata); + KapaPlotVector (kapa, nF, xFaint->data.F32); + KapaPlotVector (kapa, nF, yFaint->data.F32); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 2; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (kapa, nB, &graphdata); + KapaPlotVector (kapa, nB, xBright->data.F32); + KapaPlotVector (kapa, nB, yBright->data.F32); + + psLogMsg ("psphot", 3, "saving plot to %s", file->filename); + KapaPNG (kapa, file->filename); + + psFree (xBright); + psFree (yBright); + psFree (xFaint); + psFree (yFaint); + + return true; +} + +// plot the sx, sy, sxy as vector field, +// plot the PSF measured sx, sy, sxy as vector field +// pull the sources from the config / file? +bool psphotPlotPSFModel (pmConfig *config, pmFPAview *view, psArray *sources) { + + // select model pixels (from output background model file, or create internal file) + pmFPAfile *file = psMetadataLookupPtr (NULL, config->files, "PSPHOT.PSFMODEL.PLT"); + if (file == NULL) { + psLogMsg ("psphot", 3, "skipping psf model plot"); + return false; + } + + // pmFPAfileOpen defers disk I/O for KAPA files: just get the correct name + pmFPAfileOpen (file, view, config); + + Graphdata graphdata; + + psLogMsg ("psphot", 3, "creating psf model plot"); + + int kapa = psphotKapaOpen (); + if (kapa == -1) { + psError(PSPHOT_ERR_UNKNOWN, true, "failure to open kapa"); + return false; + } + + // XXX make the aspect-ratio match the image + KapaResize (kapa, 800, 800); + KapaInitGraph (&graphdata); + + psVector *xMNT = psVectorAllocEmpty (2*sources->n, PS_TYPE_F32); + psVector *yMNT = psVectorAllocEmpty (2*sources->n, PS_TYPE_F32); + psVector *xPSF = psVectorAllocEmpty (2*sources->n, PS_TYPE_F32); + psVector *yPSF = psVectorAllocEmpty (2*sources->n, PS_TYPE_F32); + psVector *xMIN = psVectorAllocEmpty (2*sources->n, PS_TYPE_F32); + psVector *yMIN = psVectorAllocEmpty (2*sources->n, PS_TYPE_F32); + + // construct the plot vectors + int nMNT = 0; + int nPSF = 0; + int nMIN = 0; + float dx = 0; + float dy = 0; + float scale = 10; + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (source->moments == NULL) continue; + if (source->moments->SN < 25) continue; + if (source->type != PM_SOURCE_TYPE_STAR) continue; + + pmModel *model = source->modelPSF; + if (model == NULL) continue; + + psF32 *PAR = model->params->data.F32; + + psEllipseMoments moments; + moments.x2 = source->moments->Sx; + moments.y2 = source->moments->Sy; + moments.xy = source->moments->Sxy; + + psEllipseShape shape; + shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0); + shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0); + shape.sxy = PAR[PM_PAR_SXY]; + + // force the axis ratio to be < 20.0 + psEllipseAxes axes_mnt = psEllipseMomentsToAxes (moments, 20.0); + psEllipseAxes axes_psf = psEllipseShapeToAxes (shape, 20.0); + + // moments major axis + dx = scale*axes_mnt.major*cos(axes_mnt.theta); + dy = scale*axes_mnt.major*sin(axes_mnt.theta); + xMNT->data.F32[nMNT] = PAR[PM_PAR_XPOS] - dx; + yMNT->data.F32[nMNT] = PAR[PM_PAR_YPOS] - dy; + nMNT++; + xMNT->data.F32[nMNT] = PAR[PM_PAR_XPOS] + dx; + yMNT->data.F32[nMNT] = PAR[PM_PAR_YPOS] + dy; + nMNT++; + + // psf major axis + dx = scale*axes_psf.major*cos(axes_psf.theta); + dy = scale*axes_psf.major*sin(axes_psf.theta); + xPSF->data.F32[nPSF] = PAR[PM_PAR_XPOS] - dx; + yPSF->data.F32[nPSF] = PAR[PM_PAR_YPOS] - dy; + nPSF++; + xPSF->data.F32[nPSF] = PAR[PM_PAR_XPOS] + dx; + yPSF->data.F32[nPSF] = PAR[PM_PAR_YPOS] + dy; + nPSF++; + + // minor axis (to show size) + dy = +scale*axes_psf.minor*cos(axes_psf.theta); + dx = -scale*axes_psf.minor*sin(axes_psf.theta); + xMIN->data.F32[nMIN] = PAR[PM_PAR_XPOS] - dx; + yMIN->data.F32[nMIN] = PAR[PM_PAR_YPOS] - dy; + nMIN++; + xMIN->data.F32[nMIN] = PAR[PM_PAR_XPOS] + dx; + yMIN->data.F32[nMIN] = PAR[PM_PAR_YPOS] + dy; + nMIN++; + + graphdata.xmin = PS_MIN(graphdata.xmin, PAR[PM_PAR_XPOS]); + graphdata.xmax = PS_MAX(graphdata.xmax, PAR[PM_PAR_XPOS]); + graphdata.ymin = PS_MIN(graphdata.ymin, PAR[PM_PAR_YPOS]); + graphdata.ymax = PS_MAX(graphdata.ymax, PAR[PM_PAR_YPOS]); + } + xMNT->n = yMNT->n = nMNT; + xPSF->n = yPSF->n = nPSF; + xMIN->n = yMIN->n = nMIN; + + float range; + range = graphdata.xmax - graphdata.xmin; + graphdata.xmax += 0.05*range; + graphdata.xmin -= 0.05*range; + range = graphdata.ymax - graphdata.ymin; + graphdata.ymax += 0.05*range; + graphdata.ymin -= 0.05*range; + + // XXX set the plot range to match the image + KapaSetLimits (kapa, &graphdata); + + KapaSetFont (kapa, "helvetica", 14); + KapaBox (kapa, &graphdata); + KapaSendLabel (kapa, "x (pixels)", KAPA_LABEL_XM); + KapaSendLabel (kapa, "y (pixels)", KAPA_LABEL_YM); + KapaSendLabel (kapa, "vector is major axis (scale by 20) : black are moments, blue are psf model, red is psf minor axis", KAPA_LABEL_XP); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 100; + graphdata.size = 0.3; + graphdata.style = 2; + KapaPrepPlot (kapa, nMNT, &graphdata); + KapaPlotVector (kapa, nMNT, xMNT->data.F32); + KapaPlotVector (kapa, nMNT, yMNT->data.F32); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 100; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (kapa, nPSF, &graphdata); + KapaPlotVector (kapa, nPSF, xPSF->data.F32); + KapaPlotVector (kapa, nPSF, yPSF->data.F32); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 100; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (kapa, nMIN, &graphdata); + KapaPlotVector (kapa, nMIN, xMIN->data.F32); + KapaPlotVector (kapa, nMIN, yMIN->data.F32); + + psLogMsg ("psphot", 3, "saving plot to %s", file->filename); + KapaPNG (kapa, file->filename); + + psFree (xMNT); + psFree (yMNT); + psFree (xPSF); + psFree (yPSF); + psFree (xMIN); + psFree (yMIN); + + return true; +} + +# else + +bool psphotPlotMoments (pmConfig *config, pmFPAview *view, psArray *sources) { + psLogMsg ("psphot", 3, "skipping moments plot"); + return true; +} + +bool psphotPlotPSFModel (pmConfig *config, pmFPAview *view, psArray *sources) { + psLogMsg ("psphot", 3, "skipping psf model plot"); + return true; +} + +# endif