Changeset 23344

Timestamp:

Mar 17, 2009, 11:03:01 AM (17 years ago)

Author:

beaumont

Message:

minor updates to relastroVisual

Location:

branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro

Files:

• include/relastro.h (modified) (1 diff)
• src/ConfigInit.c (modified) (1 diff)
• src/FitChip.c (modified) (5 diffs)
• src/ImageOps.c (modified) (1 diff)
• src/bcatalog.c (modified) (7 diffs)
• src/relastroVisual.c (modified) (10 diffs)
• src/relastroVisual.h (modified) (1 diff)

branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/include/relastro.h

@@ -87 +87 @@
 
 double SIGMA_LIM;
-int MIN_SRC_FIT;
+int SRC_MEAS_TOOFEW; //catalog objects wich fewer detections then this are ignored
 double MIN_ERROR;
 

branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/ConfigInit.c

@@ -20 +20 @@
 
   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

branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/FitChip.c

@@ -6 +6 @@
 # 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 +29 @@
     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 +42 @@
     dsort (values, Nscatter);
     dRmax = FIT_CHIP_NSIGMA*values[(int)(0.40*Nscatter)];
+    relastroVisualPlotScatter(values, dRmax, Nscatter);
     relastroVisualPlotRawRef(raw, ref, dRmax, Nmatch);
 
@@ -52 +54 @@
     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 +99 @@
       XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, coords);
     }
-
-    relastroVisualPlotRawRef(fit, coords, dRmax, Nmatch);
 
   }

branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/ImageOps.c

@@ -293 +293 @@
     // 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) {

branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/bcatalog.c

@@ -2 +2 @@
 
 int bcatalog (Catalog *subcatalog, Catalog *catalog) {
-  
+
   int i, j, k, offset, found;
   int NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm, Nsecfilt;
@@ -23 +23 @@
   /* 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 +47 @@
       /* 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 +75 @@
       /* 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 +89 @@
       /* 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 +110 @@
       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 +139 @@
 
   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);

branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastroVisual.c

@@ -6 +6 @@
 #include "relastro.h"
 
+#define TESTING
+
 #define KAPAX 700
 #define KAPAY 700
@@ -11 +13 @@
 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 +30 @@
       return 0;
     }
-    KapaResize (*kapid, KAPAX, KAPAY);
+    //    KapaResize (*kapid, KAPAX, KAPAY);
   }
   return 1;
@@ -46 +52 @@
 
 /** 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, &section);
+    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, &section);
+    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, &section);
+    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, &section);
+    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 +157 @@
     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 +185 @@
     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 +197 @@
     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, &section);
+
+    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 +307 @@
   }
 
-  ALLOCATE(xVec, float, numFit);
-  ALLOCATE(yVec, float, numFit);
+  if (numFit == 0) return 0;
 
   Graphdata graphdata;
@@ -184 +317 @@
   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 +354 @@
   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);

branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastroVisual.h

@@ -3 +3 @@
 
 int relastroVisualPlotRawRef(StarData *raw, StarData *ref, double dRmax, int numObj);
+int relastroVisualPlotScatter(double values[], double thresh, int npts);
 
 #endif