Index: trunk/psModules/src/objects/Makefile.am =================================================================== --- trunk/psModules/src/objects/Makefile.am (revision 20522) +++ trunk/psModules/src/objects/Makefile.am (revision 20582) @@ -39,4 +39,5 @@ pmSourcePlotMoments.c \ pmSourcePlotApResid.c \ + pmSourceVisual.c \ pmResiduals.c \ pmPSF.c \ @@ -74,4 +75,5 @@ pmSourceIO.h \ pmSourcePlots.h \ + pmSourceVisual.h \ pmResiduals.h \ pmPSF.h \ Index: trunk/psModules/src/objects/pmPSFtry.c =================================================================== --- trunk/psModules/src/objects/pmPSFtry.c (revision 20522) +++ trunk/psModules/src/objects/pmPSFtry.c (revision 20582) @@ -5,6 +5,6 @@ * @author EAM, IfA * - * @version $Revision: 1.65 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-10-29 01:04:27 $ + * @version $Revision: 1.66 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-11-08 01:52:34 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -34,4 +34,5 @@ #include "pmSourceFitModel.h" #include "pmSourcePhotometry.h" +#include "pmSourceVisual.h" bool printTrendMap (pmTrend2D *trend) { @@ -683,4 +684,5 @@ stdev = psStatsGetValue (trend->stats, stdevOption); psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0->n); + pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask); trend = psf->params->data[PM_PAR_E1]; @@ -689,4 +691,5 @@ stdev = psStatsGetValue (trend->stats, stdevOption); psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1->n); + pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask); trend = psf->params->data[PM_PAR_E2]; @@ -695,4 +698,5 @@ stdev = psStatsGetValue (trend->stats, stdevOption); psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2->n); + pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask); if (!status) { @@ -853,4 +857,5 @@ psImageMapCleanup (trend->map); // printTrendMap (trend); + pmSourceVisualPSFModelResid (trend, x, y, e0obs, mask); trend = psf->params->data[PM_PAR_E1]; @@ -862,4 +867,5 @@ psImageMapCleanup (trend->map); // printTrendMap (trend); + pmSourceVisualPSFModelResid (trend, x, y, e1obs, mask); trend = psf->params->data[PM_PAR_E2]; @@ -871,4 +877,5 @@ psImageMapCleanup (trend->map); // printTrendMap (trend); + pmSourceVisualPSFModelResid (trend, x, y, e2obs, mask); } psf->psfTrendStats->clipIter = nIter; // restore default setting Index: trunk/psModules/src/objects/pmSourceVisual.c =================================================================== --- trunk/psModules/src/objects/pmSourceVisual.c (revision 20582) +++ trunk/psModules/src/objects/pmSourceVisual.c (revision 20582) @@ -0,0 +1,256 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include "pmTrend2D.h" +#include "pmSourceVisual.h" + +# if (HAVE_KAPA) +# include + +// functions used to visualize the analysis as it goes +// these are invoked by the -visual options + +static bool isVisual = false; +static int kapa1 = -1; +// static int kapa2 = -1; +// static int kapa3 = -1; + +bool pmSourcePlotPoints3D (int myKapa, Graphdata *graphdata, psVector *xn, psVector *yn, psVector *zn, float theta, float phi); + +bool pmSourceSetVisual (bool mode) { + + isVisual = mode; + return true; +} + +bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask) { + + KapaSection section; // put the positive profile in one and the residuals in another? + + Graphdata graphdata; + + if (!isVisual) return true; + + if (kapa1 == -1) { + kapa1 = KapaOpenNamedSocket ("kapa", "pmSource:plots"); + if (kapa1 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + isVisual = false; + return false; + } + } + + KapaClearPlots (kapa1); + KapaInitGraph (&graphdata); + + float min = +1e32; + float max = -1e32; + float Min = +1e32; + float Max = -1e32; + + psVector *resid = psVectorAlloc (x->n, PS_TYPE_F32); + psVector *model = psVectorAlloc (x->n, PS_TYPE_F32); + + for (int i = 0; i < x->n; i++) { + model->data.F32[i] = pmTrend2DEval (trend, x->data.F32[i], y->data.F32[i]); + resid->data.F32[i] = param->data.F32[i] - model->data.F32[i]; + if (mask->data.U8[i]) continue; + min = PS_MIN (min, resid->data.F32[i]); + max = PS_MAX (max, resid->data.F32[i]); + Min = PS_MIN (min, param->data.F32[i]); + Max = PS_MAX (max, param->data.F32[i]); + } + + psVector *xn = psVectorAlloc (x->n, PS_TYPE_F32); + psVector *yn = psVectorAlloc (x->n, PS_TYPE_F32); + psVector *zn = psVectorAlloc (x->n, PS_TYPE_F32); + psVector *Zn = psVectorAlloc (x->n, PS_TYPE_F32); + psVector *Fn = psVectorAlloc (x->n, PS_TYPE_F32); + for (int i = 0; i < x->n; i++) { + xn->data.F32[i] = x->data.F32[i] / 5000.0; + yn->data.F32[i] = y->data.F32[i] / 5000.0; + zn->data.F32[i] = (resid->data.F32[i] - min) / (max - min); + Zn->data.F32[i] = (param->data.F32[i] - Min) / (Max - Min); + Fn->data.F32[i] = (model->data.F32[i] - Min) / (Max - Min); + } + + // view 1 on resid + section.dx = 0.5; + section.dy = 0.33; + section.x = 0.0; + section.y = 0.0; + section.name = NULL; + psStringAppend (§ion.name, "a1"); + KapaSetSection (kapa1, §ion); + psFree (section.name); + pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, zn, 30.0*PS_RAD_DEG, -15.0*PS_RAD_DEG); + + // view 2 on resid + section.dx = 0.5; + section.dy = 0.33; + section.x = 0.5; + section.y = 0.0; + section.name = NULL; + psStringAppend (§ion.name, "a2"); + KapaSetSection (kapa1, §ion); + psFree (section.name); + pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, zn, -60.0*PS_RAD_DEG, -15.0*PS_RAD_DEG); + + // view 3 on resid + section.dx = 0.5; + section.dy = 0.33; + section.x = 0.0; + section.y = 0.33; + section.name = NULL; + psStringAppend (§ion.name, "a3"); + KapaSetSection (kapa1, §ion); + psFree (section.name); + pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Zn, 30.0*PS_RAD_DEG, -15.0*PS_RAD_DEG); + + // view 4 on resid + section.dx = 0.5; + section.dy = 0.33; + section.x = 0.5; + section.y = 0.33; + section.name = NULL; + psStringAppend (§ion.name, "a4"); + KapaSetSection (kapa1, §ion); + psFree (section.name); + pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Zn, -60.0*PS_RAD_DEG, -15.0*PS_RAD_DEG); + + // view 5 on resid + section.dx = 0.5; + section.dy = 0.33; + section.x = 0.0; + section.y = 0.66; + section.name = NULL; + psStringAppend (§ion.name, "a5"); + KapaSetSection (kapa1, §ion); + psFree (section.name); + pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Fn, 30.0*PS_RAD_DEG, -15.0*PS_RAD_DEG); + + // view 6 on resid + section.dx = 0.5; + section.dy = 0.33; + section.x = 0.5; + section.y = 0.66; + section.name = NULL; + psStringAppend (§ion.name, "a6"); + KapaSetSection (kapa1, §ion); + psFree (section.name); + pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Fn, -60.0*PS_RAD_DEG, -15.0*PS_RAD_DEG); + + psFree (resid); + + psFree (xn); + psFree (yn); + psFree (zn); + psFree (Zn); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + return true; +} + +// send in normalized points +bool pmSourcePlotPoints3D (int myKapa, Graphdata *graphdata, psVector *xn, psVector *yn, psVector *zn, float theta, float phi) { + + psVector *xv = psVectorAlloc (PS_MAX(6, 2*xn->n), PS_TYPE_F32); + psVector *yv = psVectorAlloc (PS_MAX(6, 2*xn->n), PS_TYPE_F32); + psVector *zv = psVectorAlloc (PS_MAX(6, 2*xn->n), PS_TYPE_F32); + + graphdata->xmin = +1e32; + graphdata->xmax = -1e32; + graphdata->ymin = +1e32; + graphdata->ymax = -1e32; + + for (int i = 0; i < xn->n; i++) { + xv->data.F32[2*i+0] = +xn->data.F32[i]*cos(theta) + yn->data.F32[i]*sin(theta)*cos(phi) + zn->data.F32[i]*sin(theta)*sin(phi); + yv->data.F32[2*i+0] = -xn->data.F32[i]*sin(theta) + yn->data.F32[i]*cos(theta)*cos(phi) + zn->data.F32[i]*cos(theta)*sin(phi); + zv->data.F32[2*i+0] = -yn->data.F32[i]*sin(phi) + zn->data.F32[i]*cos(phi); + xv->data.F32[2*i+1] = +xn->data.F32[i]*cos(theta) + yn->data.F32[i]*sin(theta)*cos(phi); + yv->data.F32[2*i+1] = -xn->data.F32[i]*sin(theta) + yn->data.F32[i]*cos(theta)*cos(phi); + zv->data.F32[2*i+1] = -yn->data.F32[i]*sin(phi); + graphdata->xmin = PS_MIN(graphdata->xmin, xv->data.F32[2*i+0]); + graphdata->xmax = PS_MAX(graphdata->xmax, xv->data.F32[2*i+0]); + graphdata->ymin = PS_MIN(graphdata->ymin, zv->data.F32[2*i+0]); + graphdata->ymax = PS_MAX(graphdata->ymax, zv->data.F32[2*i+0]); + graphdata->xmin = PS_MIN(graphdata->xmin, xv->data.F32[2*i+1]); + graphdata->xmax = PS_MAX(graphdata->xmax, xv->data.F32[2*i+1]); + graphdata->ymin = PS_MIN(graphdata->ymin, zv->data.F32[2*i+1]); + graphdata->ymax = PS_MAX(graphdata->ymax, zv->data.F32[2*i+1]); + } + xv->n = xn->n; + + // examine sources to set data range + KapaSetLimits (myKapa, graphdata); + + // KapaSetFont (myKapa, "helvetica", 14); + // KapaBox (myKapa, graphdata); + // KapaSendLabel (myKapa, "&ss&h_x| (pixels)", KAPA_LABEL_XM); + // KapaSendLabel (myKapa, "&ss&h_y| (pixels)", KAPA_LABEL_YM); + + graphdata->color = KapaColorByName ("black"); + graphdata->ptype = 100; + graphdata->size = 0.5; + graphdata->style = 2; + KapaPrepPlot (myKapa, xv->n, graphdata); + KapaPlotVector (myKapa, xv->n, xv->data.F32, "x"); + KapaPlotVector (myKapa, xv->n, zv->data.F32, "y"); + + graphdata->color = KapaColorByName ("blue"); + graphdata->ptype = 0; + graphdata->size = 1.5; + graphdata->style = 2; + KapaPrepPlot (myKapa, xv->n, graphdata); + KapaPlotVector (myKapa, xv->n, xv->data.F32, "x"); + KapaPlotVector (myKapa, xv->n, zv->data.F32, "y"); + + xv->n = 6; + + // set the three axis lines + xv->data.F32[0] = +0.0*cos(theta) + 0.0*sin(theta)*cos(phi) + 0.0*sin(theta)*sin(phi); + yv->data.F32[0] = -0.0*sin(theta) + 0.0*cos(theta)*cos(phi) + 0.0*cos(theta)*sin(phi); + zv->data.F32[0] = - 0.0*sin(phi) + 0.0*cos(phi); + xv->data.F32[1] = +1.0*cos(theta) + 0.0*sin(theta)*cos(phi); + yv->data.F32[1] = -1.0*sin(theta) + 0.0*cos(theta)*cos(phi); + zv->data.F32[1] = - 0.0*sin(phi) + 0.0*cos(phi); + + xv->data.F32[2] = +0.0*cos(theta) + 0.0*sin(theta)*cos(phi) + 0.0*sin(theta)*sin(phi); + yv->data.F32[2] = -0.0*sin(theta) + 0.0*cos(theta)*cos(phi) + 0.0*cos(theta)*sin(phi); + zv->data.F32[2] = - 0.0*sin(phi) + 0.0*cos(phi); + xv->data.F32[3] = +0.0*cos(theta) + 1.0*sin(theta)*cos(phi); + yv->data.F32[3] = -0.0*sin(theta) + 1.0*cos(theta)*cos(phi); + zv->data.F32[3] = - 1.0*sin(phi) + 0.0*cos(phi); + + xv->data.F32[4] = +0.0*cos(theta) + 0.0*sin(theta)*cos(phi) + 1.0*sin(theta)*sin(phi); + yv->data.F32[4] = -0.0*sin(theta) + 0.0*cos(theta)*cos(phi) + 1.0*cos(theta)*sin(phi); + zv->data.F32[4] = - 0.0*sin(phi) + 1.0*cos(phi); + xv->data.F32[5] = +0.0*cos(theta) + 0.0*sin(theta)*cos(phi); + yv->data.F32[5] = -0.0*sin(theta) + 0.0*cos(theta)*cos(phi); + zv->data.F32[5] = - 0.0*sin(phi) + 0.0*cos(phi); + + graphdata->color = KapaColorByName ("red"); + graphdata->ptype = 100; + graphdata->size = 0.5; + graphdata->style = 2; + KapaPrepPlot (myKapa, xv->n, graphdata); + KapaPlotVector (myKapa, xv->n, xv->data.F32, "x"); + KapaPlotVector (myKapa, xv->n, zv->data.F32, "y"); + + psFree (xv); + psFree (yv); + psFree (zv); + + return true; +} + +# else +# endif Index: trunk/psModules/src/objects/pmSourceVisual.h =================================================================== --- trunk/psModules/src/objects/pmSourceVisual.h (revision 20582) +++ trunk/psModules/src/objects/pmSourceVisual.h (revision 20582) @@ -0,0 +1,21 @@ +/* @file pmKapaPlots.h + * @brief functions to make plots with the external program 'kapa' + * + * @author EAM, IfA + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-11-08 01:52:34 $ + * Copyright 2006 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_SOURCE_VISUAL_H +#define PM_SOURCE_VISUAL_H + +/// @addtogroup Extras Miscellaneous Funtions +/// @{ + +bool pmSourceSetVisual (bool mode); +bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask); + +/// @} +#endif // PM_KAPA_PLOTS_H