Index: /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/include/relastro.h =================================================================== --- /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/include/relastro.h (revision 23343) +++ /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/include/relastro.h (revision 23344) @@ -87,5 +87,5 @@ double SIGMA_LIM; -int MIN_SRC_FIT; +int SRC_MEAS_TOOFEW; //catalog objects wich fewer detections then this are ignored double MIN_ERROR; Index: /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/ConfigInit.c (revision 23343) +++ /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/ConfigInit.c (revision 23344) @@ -20,5 +20,5 @@ GetConfig (config, "RELASTRO_SIGMA_LIM", "%lf", 0, &SIGMA_LIM); - GetConfig (config, "RELASTRO_MIN_SRC_FIT", "%d", 0, &MIN_SRC_FIT); + GetConfig (config, "RELASTRO_SRC_MEAS_TOOFEW", "%d", 0, &SRC_MEAS_TOOFEW); // XXX these are used in relphot to identify poor stars / images -- define Index: /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/FitChip.c =================================================================== --- /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/FitChip.c (revision 23343) +++ /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/FitChip.c (revision 23344) @@ -6,4 +6,5 @@ # define FIT_CHIP_NITER 3 # define FIT_CHIP_NSIGMA 3.0 + // XXX we should test if the fit is sufficiently constrained across the chip, or if the @@ -28,5 +29,5 @@ for (i = Nscatter = 0; i < Nmatch; i++) { if (raw[i].mask) continue; - if (raw[i].dMag > FIT_CHIP_MAX_ERROR) continue; + if (isnan(raw[i].dMag) || raw[i].dMag > FIT_CHIP_MAX_ERROR) continue; dL = raw[i].L - ref[i].L; @@ -41,4 +42,5 @@ dsort (values, Nscatter); dRmax = FIT_CHIP_NSIGMA*values[(int)(0.40*Nscatter)]; + relastroVisualPlotScatter(values, dRmax, Nscatter); relastroVisualPlotRawRef(raw, ref, dRmax, Nmatch); @@ -52,4 +54,5 @@ for (i = 0; i < Nmatch; i++) { if (raw[i].mask) continue; + if (isnan(raw[i].dMag) || raw[i].dMag > FIT_CHIP_MAX_ERROR) continue; // only keep objects within dRmax @@ -96,6 +99,4 @@ XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, coords); } - - relastroVisualPlotRawRef(fit, coords, dRmax, Nmatch); } Index: /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/ImageOps.c =================================================================== --- /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/ImageOps.c (revision 23343) +++ /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/ImageOps.c (revision 23344) @@ -293,7 +293,10 @@ // an object with only one detection provides no information about the image calibration raw[i].mask = FALSE; - if (catalog[c].average[n].Nmeasure < MIN_SRC_FIT) { - raw[i].mask = TRUE; - } + int mask = FALSE; + if (catalog[c].average[n].Nmeasure <= SRC_MEAS_TOOFEW) { + mask = TRUE; + } + raw[i].mask = mask; + switch (mode) { Index: /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/bcatalog.c =================================================================== --- /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/bcatalog.c (revision 23343) +++ /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/bcatalog.c (revision 23344) @@ -2,5 +2,5 @@ int bcatalog (Catalog *subcatalog, Catalog *catalog) { - + int i, j, k, offset, found; int NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm, Nsecfilt; @@ -23,5 +23,5 @@ /* exclude stars not in range or with too few measurements */ for (i = 0; i < catalog[0].Naverage; i++) { - if (catalog[0].average[i].Nmeasure < 2) continue; + if (catalog[0].average[i].Nmeasure <= SRC_MEAS_TOOFEW) continue; /* start with all stars good */ @@ -47,24 +47,24 @@ /* select measurements by photcode, or equiv photcode, if specified */ if (NphotcodesKeep > 0) { - found = FALSE; - for (k = 0; (k < NphotcodesKeep) && !found; k++) { - if (photcodesKeep[k][0].code == catalog[0].measure[offset].photcode) found = TRUE; - if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE; - } - if (!found) continue; + found = FALSE; + for (k = 0; (k < NphotcodesKeep) && !found; k++) { + if (photcodesKeep[k][0].code == catalog[0].measure[offset].photcode) found = TRUE; + if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE; + } + if (!found) continue; } if (NphotcodesSkip > 0) { - found = FALSE; - for (k = 0; (k < NphotcodesSkip) && !found; k++) { - if (photcodesSkip[k][0].code == catalog[0].measure[offset].photcode) found = TRUE; - if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE; - } - if (found) continue; + found = FALSE; + for (k = 0; (k < NphotcodesSkip) && !found; k++) { + if (photcodesSkip[k][0].code == catalog[0].measure[offset].photcode) found = TRUE; + if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE; + } + if (found) continue; } /* select measurements by time */ if (TimeSelect) { - if (catalog[0].measure[offset].t < TSTART) continue; - if (catalog[0].measure[offset].t > TSTOP) continue; + if (catalog[0].measure[offset].t < TSTART) continue; + if (catalog[0].measure[offset].t > TSTOP) continue; } @@ -75,11 +75,11 @@ /* select measurements by quality */ if (PhotFlagSelect) { - if (PhotFlagBad) { - mask = PhotFlagBad; - } else { - code = GetPhotcodebyCode (catalog[0].measure[offset].photcode); - mask = code[0].astromBadMask; - } - if (mask & catalog[0].measure[offset].photFlags) continue; + if (PhotFlagBad) { + mask = PhotFlagBad; + } else { + code = GetPhotcodebyCode (catalog[0].measure[offset].photcode); + mask = code[0].astromBadMask; + } + if (mask & catalog[0].measure[offset].photFlags) continue; } @@ -89,7 +89,7 @@ /* select measurements by mag limit */ if (ImagSelect) { - mag = PhotInst (&catalog[0].measure[offset]); - if (mag < ImagMin) continue; - if (mag > ImagMax) continue; + mag = PhotInst (&catalog[0].measure[offset]); + if (mag < ImagMin) continue; + if (mag > ImagMax) continue; } @@ -110,13 +110,13 @@ subcatalog[0].measure[Nmeasure] = catalog[0].measure[offset]; subcatalog[0].measure[Nmeasure].averef = Naverage; - if (RESET) { - subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM; - subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_AREA; + if (RESET) { + subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM; + subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_AREA; } Nmeasure ++; Nm ++; if (Nmeasure == NMEASURE) { - NMEASURE += 1000; - REALLOCATE (subcatalog[0].measure, Measure, NMEASURE); + NMEASURE += 1000; + REALLOCATE (subcatalog[0].measure, Measure, NMEASURE); } } @@ -139,6 +139,6 @@ if (VERBOSE) { - fprintf (stderr, "%d: using %d stars (%d measures) for catalog\n", i, - subcatalog[0].Naverage, subcatalog[0].Nmeasure); + fprintf (stderr, "%d: using %d stars (%d measures) for catalog\n", i, + subcatalog[0].Naverage, subcatalog[0].Nmeasure); } return (TRUE); Index: /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastroVisual.c =================================================================== --- /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastroVisual.c (revision 23343) +++ /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastroVisual.c (revision 23344) @@ -6,4 +6,6 @@ #include "relastro.h" +#define TESTING + #define KAPAX 700 #define KAPAY 700 @@ -11,7 +13,11 @@ static int kapa = -1; static int kapa2 = -1; - -static int isVisual = 1; -static int plotRawRef = 1; +static int kapa3 = -1; + +static int isVisual = FALSE; +static int plotRawRef = TRUE; +static int plotScatter = TRUE; +static int plotResid = TRUE; +static int plotVector = TRUE; /** Spawn a kapa window */ @@ -24,5 +30,5 @@ return 0; } - KapaResize (*kapid, KAPAX, KAPAY); + // KapaResize (*kapid, KAPAX, KAPAY); } return 1; @@ -46,29 +52,104 @@ /** Size graphdata to encompass all points */ -static void scaleGraphdata(float x[], float y[], Graphdata *graphdata, int n) { - float xlo = 9999999, xhi = -9999999, ylo = 9999999, yhi = -9999999; - int i; - for(i = 0; i < n; i++) { - if(x[i] < xlo) xlo = x[i]; - if(x[i] > xhi) xhi = x[i]; - if(y[i] < ylo) ylo = y[i]; - if(y[i] > yhi) yhi = y[i]; - } - graphdata->xmin = xlo; - graphdata->ymin = ylo; - graphdata->xmax = xhi; - graphdata->ymax = yhi; -} +static int scaleGraphdata(float x[], float y[], Graphdata *graphdata, int n) { + float xlo = FLT_MAX, xhi = -FLT_MAX, ylo = FLT_MAX, yhi = -FLT_MAX; + int i; + int goodData = 0; + for(i = 0; i < n; i++) { + goodData++; + if(!finite(x[i]) || !finite(y[i])) continue; + if(x[i] < xlo) xlo = x[i]; + if(x[i] > xhi) xhi = x[i]; + if(y[i] < ylo) ylo = y[i]; + if(y[i] > yhi) yhi = y[i]; + } + if (goodData < 2) return 0; + graphdata->xmin = xlo; + graphdata->ymin = ylo; + graphdata->xmax = xhi; + graphdata->ymax = yhi; + return 1; +} + +static int residPlot(float x[], float y[], + float xVec[], float yVec[], + int npts, int *kapaID) { + if (!plotResid) return 1; + Graphdata graphdata; + KapaSection section; + + if (!initWindow(kapaID)) return 0; + + KapaInitGraph(&graphdata); + KapaClearPlots(*kapaID); + KapaSetFont(*kapaID, "helvetica", 14); + + section.name = "0"; + section.x = 0.0; section.y = 0.0; + section.dx = .45, section.dy = .45; + graphdata.ptype = 7; + graphdata.style = 2; + + KapaSetSection(*kapaID, §ion); + if(!scaleGraphdata(x, xVec, &graphdata, npts)) return 0; + KapaSetLimits(*kapaID, &graphdata); + KapaBox(*kapaID, &graphdata); + KapaSendLabel(*kapaID, "x", KAPA_LABEL_XM); + KapaSendLabel(*kapaID, "dx", KAPA_LABEL_YM); + KapaPrepPlot(*kapaID, npts, &graphdata); + KapaPlotVector(*kapaID, npts, x, "x"); + KapaPlotVector(*kapaID, npts, xVec, "y"); + + section.x = .5; section.y = 0; section.name="1"; + KapaSetSection(*kapaID, §ion); + if(!scaleGraphdata(x, yVec, &graphdata, npts)) return 0; + KapaSetLimits(*kapaID, &graphdata); + KapaBox(*kapaID, &graphdata); + KapaSendLabel(*kapaID, "x", KAPA_LABEL_XM); + KapaSendLabel(*kapaID, "dy", KAPA_LABEL_YP); + KapaPrepPlot(*kapaID, npts, &graphdata); + KapaPlotVector(*kapaID, npts, x, "x"); + KapaPlotVector(*kapaID, npts, yVec, "y"); + + section.x = .0; section.y = .5; section.name="2"; + KapaSetSection(*kapaID, §ion); + if(!scaleGraphdata(y, xVec, &graphdata, npts)) return 0;; + KapaSetLimits(*kapaID, &graphdata); + KapaBox(*kapaID, &graphdata); + KapaSendLabel(*kapaID, "y", KAPA_LABEL_XM); + KapaSendLabel(*kapaID, "dx", KAPA_LABEL_YM); + KapaPrepPlot(*kapaID, npts, &graphdata); + KapaPlotVector(*kapaID, npts, y, "x"); + KapaPlotVector(*kapaID, npts, xVec, "y"); + + section.x = .5; section.y = .5; section.name="3"; + KapaSetSection(*kapaID, §ion); + if(!scaleGraphdata(y, yVec, &graphdata, npts)) return 0; + KapaSetLimits(*kapaID, &graphdata); + KapaBox(*kapaID, &graphdata); + KapaSendLabel(*kapaID, "y", KAPA_LABEL_XM); + KapaSendLabel(*kapaID, "dy", KAPA_LABEL_YP); + KapaPrepPlot(*kapaID, npts, &graphdata); + KapaPlotVector(*kapaID, npts, y, "x"); + KapaPlotVector(*kapaID, npts, yVec, "y"); + + return 1; +} + /**Plot a vector field*/ static int plotVectorField(float x[], float y[], float xVec[], float yVec[], - int npts, int *kapaID) { - - Graphdata graphdata; - float singleX[2], singleY[2]; - float maxVecLength, vecScaleFactor; - float graphSize; - int i; + int npts, int *kapaID, double maxVecLength) { + + if(!plotVector) return 1; + + Graphdata graphdata; + float singleX[2], singleY[2]; + float vecScaleFactor; + float graphSize; + int i; + char plotTitle[50]; + sprintf(plotTitle, "Maximum Vector Size = %5.1e", maxVecLength); if (!initWindow(kapaID)) return 0; @@ -76,26 +157,23 @@ KapaInitGraph(&graphdata); KapaClearPlots(*kapaID); - scaleGraphdata(x, y, &graphdata, npts); - - //scale the raw xvec, yvecs down to fit into the window - maxVecLength = 0; - for(i = 0; i < npts; i++) { - if (fabs(xVec[i]) > maxVecLength) { - maxVecLength = fabs(xVec[i]); - } - if (fabs(yVec[i]) > maxVecLength) { - maxVecLength = fabs(yVec[i]); - } - } + if(!scaleGraphdata(x, y, &graphdata, npts)) return 0; + graphSize = graphdata.xmax - graphdata.xmin; if ((graphdata.ymax - graphdata.ymin) > graphSize) { graphSize = graphdata.ymax - graphdata.ymin; } - vecScaleFactor = graphSize * 0.1 / maxVecLength; + + vecScaleFactor = graphSize * 0.02 / (float)maxVecLength; +#ifdef TESTING + fprintf(stderr, "GraphSize: %e\n maxVecLength: %e\n vecScaleFactor: %e\n", + graphSize, maxVecLength, vecScaleFactor); +#endif + KapaSetFont (*kapaID, "helvetica", 14); KapaSetLimits(*kapaID, &graphdata); KapaBox(*kapaID, &graphdata); - + KapaSendLabel(*kapaID, plotTitle, KAPA_LABEL_XP); + graphdata.ptype = 7; graphdata.style = 2; @@ -107,10 +185,10 @@ graphdata.ptype = 0; graphdata.style = 0; + graphdata.color = KapaColorByName("blue"); for(i = 0; i < npts; i++) { singleX[0] = x[i]; singleY[0] = y[i]; - singleX[1] = xVec[i] * vecScaleFactor; - singleY[1] = yVec[i] * vecScaleFactor; - + singleX[1] = x[i] + xVec[i] * vecScaleFactor; + singleY[1] = y[i] + yVec[i] * vecScaleFactor; KapaPrepPlot(*kapaID, 2, &graphdata); KapaPlotVector(*kapaID, 2, singleX, "x"); @@ -119,4 +197,60 @@ return 1; } + +int relastroVisualPlotScatter(double values[], double thresh, int npts) { + float *x, *data; + int i; + float xline[2], yline[2]; + if (!isVisual || !plotScatter) return 1; + if (!initWindow(&kapa2)) return 0; + + ALLOCATE(x, float, npts); + ALLOCATE(data, float, npts); + + for(i = 0; i < npts; i++) { + x[i] = i; + data[i] = (float) values[i]; + } + + Graphdata graphdata; + KapaSection section; + section.x = 0; section.y = 0; section.dx = 1; section.dy = 1; + section.name = "junk"; + + KapaInitGraph(&graphdata); + KapaClearPlots(kapa2); + KapaSetSection(kapa2, §ion); + + graphdata.ptype = 0; + graphdata.style = 0; + graphdata.xmin = 0; + graphdata.xmax = npts; + graphdata.ymin = 0; + graphdata.ymax = data[npts-1]; + + KapaSetFont(kapa2, "helvetica", 14); + KapaSetLimits(kapa2, &graphdata); + KapaBox(kapa2, &graphdata); + KapaSendLabel( kapa2, "Object", KAPA_LABEL_XM); + KapaSendLabel( kapa2, "Offset(pixels)", KAPA_LABEL_YM); + KapaSendLabel( kapa2, "Astrometric Offset with fit cutoff", + KAPA_LABEL_XP); + KapaPrepPlot(kapa2, npts, &graphdata); + KapaPlotVector(kapa2, npts, x, "x"); + KapaPlotVector(kapa2, npts, data, "y"); + + graphdata.color = KapaColorByName("red"); + KapaPrepPlot(kapa2, 2, &graphdata); + yline[0] = (float) thresh; + yline[1] = (float) thresh; + xline[0] = graphdata.xmin; + xline[1] = graphdata.xmax; + KapaPlotVector(kapa2, 2, xline, "x"); + KapaPlotVector(kapa2, 2, yline, "y"); + + askUser(&plotScatter); + return 1; +} + @@ -173,6 +307,5 @@ } - ALLOCATE(xVec, float, numFit); - ALLOCATE(yVec, float, numFit); + if (numFit == 0) return 0; Graphdata graphdata; @@ -184,14 +317,8 @@ graphdata.style = 2; - scaleGraphdata(rawXfit, rawYfit, &graphdata, numFit); - // scaleGraphdata(refXfit, refYfit, &graphdata, numFit); - - //plot the fitted objects as vectors - for(i = 0; i < numFit; i++) { - xVec[i] = rawXfit[i] - refXfit[i]; - yVec[i] = rawYfit[i] - refYfit[i]; - } - - plotVectorField(refX, refY, xVec, yVec, numFit, &kapa2); + if(!scaleGraphdata(rawXfit, rawYfit, &graphdata, numFit)) { + fprintf(stderr, "Not enough finite points for plotting"); + return 0; + } KapaSetFont(kapa, "helvetica", 14); @@ -227,4 +354,21 @@ KapaPlotVector(kapa, numFit, refYfit, "y"); + ALLOCATE(xVec, float, numFit); + ALLOCATE(yVec, float, numFit); + + //plot the fitted objects as vectors + for(i = 0; i < numFit; i++) { + xVec[i] = rawXfit[i] - refXfit[i]; + yVec[i] = rawYfit[i] - refYfit[i]; + } + + plotVectorField(rawXfit, rawYfit, xVec, yVec, numFit, &kapa3, dRmax); + if(!residPlot(rawXfit, rawYfit, xVec, yVec, numFit, &kapa2)) { + fprintf(stderr, "Unable to plot residuals"); + return 0; + } + + FREE(xVec); + FREE(yVec); askUser(&plotRawRef); Index: /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastroVisual.h =================================================================== --- /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastroVisual.h (revision 23343) +++ /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastroVisual.h (revision 23344) @@ -3,4 +3,5 @@ int relastroVisualPlotRawRef(StarData *raw, StarData *ref, double dRmax, int numObj); +int relastroVisualPlotScatter(double values[], double thresh, int npts); #endif