Changeset 30616 for trunk/Ohana/src/relastro/src/UpdateChips.c

Timestamp:

Feb 13, 2011, 11:30:02 AM (15 years ago)

Author:

eugene

Message:

various fixes to photdbc, relastro and relphot based on analysis of the first reference database

File:

• trunk/Ohana/src/relastro/src/UpdateChips.c (modified) (3 diffs)

trunk/Ohana/src/relastro/src/UpdateChips.c

@@ -1 +1 @@
 # include "relastro.h"
+int plotChipFits (double *Ro, double *Do, char *mode, int Nimage);
+int saveCenter (Image *image, double *Ro, double *Do, int im);
 
 int UpdateChips (Catalog *catalog, int Ncatalog) {
 
+  int Nskip, Nmosaic, NnewFit, NoldFit;
+
   /* we can measure new image parameters for each non-mosaic chip independently */
-  off_t i, Nimage, Nraw, Nref;
+  off_t i, Nimage, Nraw, Nref, nFitAstr;
   Image *image;
   StarData *raw, *ref;
   Coords *oldCoords;
+  float dXpixSys, dYpixSys;
+  double *Ro, *Do;
+  char *mode;
+
+  Nskip = Nmosaic = NnewFit = NoldFit = 0;
 
   image = getimages (&Nimage);
 
+  // save fit results for summary plot
+  ALLOCATE (Ro, double, Nimage);
+  ALLOCATE (Do, double, Nimage);
+  ALLOCATE (mode, char, Nimage);
+
   for (i = 0; i < Nimage; i++) {
 
     /* skip all except WRP images */
-    if (strcmp(&image[i].coords.ctype[4], "-WRP")) continue;
+    if (strcmp(&image[i].coords.ctype[4], "-WRP")) {
+      Nmosaic ++;
+      mode[i] = 0;
+      continue;
+    }
 
     /* convert measure coordinates to raw entries */
     raw = getImageRaw (catalog, Ncatalog, i, &Nraw, MODE_MOSAIC);
-    if (!raw) continue;
+    if (!raw) {
+      Nskip ++;
+      mode[i] = 0;
+      continue;
+    }
 
     /* convert average coordinates to ref entries */
     ref = getImageRef (catalog, Ncatalog, i, &Nref, MODE_MOSAIC);
-    if (!ref) continue;
+    if (!ref) {
+      Nskip ++;
+      mode[i] = 0;
+      continue;
+    }
 
     // note that Nraw & Nref must be equal: if not, we made a programming error in one of these two functions.
     assert (Nraw == Nref);
 
+    // save these in case of failure
     saveCoords (&image[i].coords, i);
+    dXpixSys = image[i].dXpixSys;
+    dYpixSys = image[i].dYpixSys;
+    nFitAstr = image[i].nFitAstrom;
 
     // FitChip does iterative, clipped fitting
@@ -33 +63 @@
     if (!FitChip (raw, ref, Nraw, &image[i])) {
       if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name,  i,  Nraw);
+
+      // restore status quo ante
       oldCoords = getCoords (i);
       memcpy (&image[i].coords, oldCoords, sizeof(Coords));
+      image[i].dXpixSys = dXpixSys;
+      image[i].dYpixSys = dYpixSys; 
+      image[i].nFitAstrom = nFitAstr;
+
+      saveCenter (image, &Ro[i], &Do[i], i);
+      mode[i] = 1;
+      NoldFit ++;
       free (raw);
       free (ref);
@@ -42 +81 @@
     if (!checkStarMap (i)) {
       if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name,  i,  Nraw);
+      // restore status quo ante
       oldCoords = getCoords (i);
       memcpy (&image[i].coords, oldCoords, sizeof(Coords));
+      image[i].dXpixSys = dXpixSys;
+      image[i].dYpixSys = dYpixSys; 
+      image[i].nFitAstrom = nFitAstr;
+
+      saveCenter (image, &Ro[i], &Do[i], i);
+      mode[i] = 2;
       image[i].flags |= ID_IMAGE_ASTROM_POOR;
-    }
-
+      NoldFit ++;
+      free (raw);
+      free (ref);
+      continue;
+    } 
+
+    saveCenter (image, &Ro[i], &Do[i], i);
+    mode[i] = 3;
+    NnewFit ++;
     free (raw);
     free (ref);
   }
 
+  plotChipFits (Ro, Do, mode, Nimage);
+
+  fprintf (stderr, "UpdateChips: %d fitted, %d keep old, %d skipped, %d mosaic (skipped)\n", NnewFit, NoldFit, Nskip, Nmosaic);
   return (TRUE);
 }
 
+int saveCenter (Image *image, double *Ro, double *Do, int im) {
+
+  Mosaic *mosaic;
+  Coords *moscoords, *imcoords;
+  double X, Y, L, M, P, Q, R, D;
+
+  moscoords = NULL;
+  if (!strcmp(&image[im].coords.ctype[4], "-WRP")) {
+    mosaic = getMosaicForImage (im);
+    if (mosaic == NULL) return FALSE;  // if we cannot find the associated image, skip it
+    moscoords = &mosaic[0].coords;
+  }
+  imcoords = &image[im].coords;
+  
+  if (!strcmp(&image[im].coords.ctype[4], "-WRP")) {
+    X = 0.5*image[im].NX; 
+    Y = 0.5*image[im].NY;
+  } else {
+    X = 0.0; 
+    Y = 0.0;
+  }
+
+  if (moscoords == NULL) {
+    // this is a Simple image (not a mosaic)
+    // note that for a Simple image, L,M = P,Q
+    XY_to_LM (&L, &M, X, Y, imcoords);
+    LM_to_RD (&R, &D, L, M, imcoords);
+  } else {
+    XY_to_LM (&L, &M, X, Y, imcoords);
+    XY_to_LM (&P, &Q, L, M, moscoords);
+    LM_to_RD (&R, &D, P, Q, moscoords);
+  }
+
+  double Rmid;
+  if (UserCatalog) {
+      Rmid = UserCatalogRA;
+  } else {
+      Rmid = 0.5*(UserPatch.Rmin + UserPatch.Rmax);
+  }
+
+  R = ohana_normalize_angle_to_midpoint (R, Rmid);
+
+  *Ro = R;
+  *Do = D;
+
+  return (TRUE);
+}
+ 
+int plotChipFits (double *Ro, double *Do, char *mode, int Nimage) {
+
+  static int kapa = -1; 
+
+  int i, N;
+  double Rmin, Rmax, Dmin, Dmax;
+  float *xvec, *yvec;
+  Graphdata graphdata;
+
+  if (!relastroGetVisual()) return (TRUE);
+
+  if (kapa == -1) {
+    kapa = KapaOpenNamedSocket("kapa", "relastro");
+    if (kapa == -1) {
+      fprintf (stderr, "can't open kapa window\n");
+      return FALSE;
+    }
+  }
+
+  Rmin = +720;
+  Rmax = -720;
+  Dmin = +90;
+  Dmax = -90;
+
+  // find the R, D range
+  for (i = 0; i < Nimage; i++) {
+    if (!mode[i]) continue;
+
+    Rmin = MIN(Rmin, Ro[i]);
+    Rmax = MAX(Rmax, Ro[i]);
+    Dmin = MIN(Dmin, Do[i]);
+    Dmax = MAX(Dmax, Do[i]);
+  }
+
+  ALLOCATE (xvec, float, Nimage);
+  ALLOCATE (yvec, float, Nimage);
+
+  bzero (&graphdata, sizeof(Graphdata));
+  plot_defaults (&graphdata);
+  graphdata.xmin = Rmin;
+  graphdata.xmax = Rmax;
+  graphdata.ymin = Dmin;
+  graphdata.ymax = Dmax;
+  graphdata.style = 2;
+  graphdata.size = 1;
+
+  KapaSetFont (kapa, "helvetica", 14);
+  KapaSetLimits (kapa, &graphdata);
+  KapaBox (kapa, &graphdata);
+
+  // *** good images ***
+  N = 0;
+  for (i = 0; i < Nimage; i++) {
+    if (mode[i] != 3) continue;
+    xvec[N] = Ro[i];
+    yvec[N] = Do[i];
+    N++;
+  }
+  graphdata.ptype = 7;
+  graphdata.color = KapaColorByName("black");
+  KapaPrepPlot (kapa, N, &graphdata);
+  KapaPlotVector (kapa, N, xvec, "x");
+  KapaPlotVector (kapa, N, yvec, "y");
+  
+  // *** reject fit ***
+  N = 0;
+  for (i = 0; i < Nimage; i++) {
+    if (mode[i] != 1) continue;
+    xvec[N] = Ro[i];
+    yvec[N] = Do[i];
+    N++;
+  }
+  graphdata.ptype = 3;
+  graphdata.color = KapaColorByName("red");
+  KapaPrepPlot (kapa, N, &graphdata);
+  KapaPlotVector (kapa, N, xvec, "x");
+  KapaPlotVector (kapa, N, yvec, "y");
+  
+  // *** divergent fit ***
+  N = 0;
+  for (i = 0; i < Nimage; i++) {
+    if (mode[i] != 2) continue;
+    xvec[N] = Ro[i];
+    yvec[N] = Do[i];
+    N++;
+  }
+  graphdata.ptype = 2;
+  graphdata.color = KapaColorByName("blue");
+  KapaPrepPlot (kapa, N, &graphdata);
+  KapaPlotVector (kapa, N, xvec, "x");
+  KapaPlotVector (kapa, N, yvec, "y");
+  
+  free (xvec);
+  free (yvec);
+
+  return (TRUE);
+}
+
+// XXX if (!FindMosaicForImage (image, Nimage, i)) { }