Changeset 37531

Timestamp:

Oct 27, 2014, 2:44:02 PM (12 years ago)

Author:

eugene

Message:

do NOT normalize RA in XY_to_RD, deal with normalization outside; fix various fakeastro issues at 0,360

Location:

branches/eam_branches/ipp-20140904/Ohana/src

Files:

• dvomerge/src/dvorepairCPT.c (modified) (1 diff)
• dvomerge/src/dvorepairDeleteImageList.c (modified) (1 diff)
• dvomerge/src/dvorepairFixTables.c (modified) (1 diff)
• fakeastro/src/fit_fake_stars.c (modified) (3 diffs)
• fakeastro/src/get_image_patch.c (modified) (3 diffs)
• fakeastro/src/make_fake_images.c (modified) (1 diff)
• fakeastro/src/make_fake_stars.c (modified) (2 diffs)
• fakeastro/src/make_fake_stars_catalog.c (modified) (3 diffs)
• libdvo/src/BoundaryTree.c (modified) (1 diff)
• libdvo/src/coordops.c (modified) (2 diffs)
• libdvo/src/skyregion_ops.c (modified) (2 diffs)
• opihi/dvo/objectcoverage.c (modified) (1 diff)
• opihi/dvo/skycoverage.c (modified) (1 diff)
• relastro/src/UpdateObjects.c (modified) (2 diffs)
• relphot/src/MosaicOps.c (modified) (1 diff)

branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/src/dvorepairCPT.c

@@ -119 +119 @@
     thisImage = MatchImage (image, Nimage, measure[i].t, measure[i].photcode, measure[i].imageID);
     XY_to_RD (&average[Nave].R, &average[Nave].D, measure[i].Xccd, measure[i].Yccd, &thisImage[0].coords);
+    average[Nave].R = ohana_normalize_angle (average[Nave].R);
 
     average[Nave].Nmeasure = 1;

branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/src/dvorepairDeleteImageList.c

@@ -400 +400 @@
     thisImage = MatchImage (image, Nimage, measure[i].t, measure[i].photcode, measure[i].imageID);
     XY_to_RD (&average[Nave].R, &average[Nave].D, measure[i].Xccd, measure[i].Yccd, &thisImage[0].coords);
+    average[Nave].R = ohana_normalize_angle (average[Nave].R);
 
     average[Nave].Nmeasure = 1;

branches/eam_branches/ipp-20140904/Ohana/src/dvomerge/src/dvorepairFixTables.c

@@ -259 +259 @@
     thisImage = MatchImage (image, Nimage, measure[i].t, measure[i].photcode, measure[i].imageID);
     XY_to_RD (&average[Nave].R, &average[Nave].D, measure[i].Xccd, measure[i].Yccd, &thisImage[0].coords);
+    average[Nave].R = ohana_normalize_angle (average[Nave].R);
 
     average[Nave].Nmeasure = 1;

branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fit_fake_stars.c

@@ -25 +25 @@
   if (Nstars >  24) order_use = 2; //  6 stars per cell
   if (Nstars >  63) order_use = 3; //  7 stars per cell
-  if (Nstars > 128) order_use = 4; //  8 stars per cell
-  if (Nstars > 225) order_use = 5; //  9 stars per cell
-  if (Nstars > 360) order_use = 6; // 10 stars per cell
+  // if (Nstars > 128) order_use = 4; //  8 stars per cell
+  // if (Nstars > 225) order_use = 5; //  9 stars per cell
+  // if (Nstars > 360) order_use = 6; // 10 stars per cell
   if (order_use < 1) {
     fprintf (stderr, "insufficient measurements (%d) for linear fit\n", Nstars);
@@ -86 +86 @@
     XY_to_LM (&Le, &Me, 0.0, 4850.0, &image[0].coords);
 
-    fprintf (stderr, "%s : (dL,dM)_1 : %f : (dL,dM)_2 : %f\n", image->name, hypot(Le-Lo, Me-Mo), hypot(Ls-Lo, Ms-Mo));
+    if (FALSE) fprintf (stderr, "%s : (dL,dM)_1 : %f : (dL,dM)_2 : %f\n", image->name, hypot(Le-Lo, Me-Mo), hypot(Ls-Lo, Ms-Mo));
   }
 
@@ -134 +134 @@
     // if we have fitted the map above, Npolyterms needs to be -1 here:
     XY_to_RD (&stars[i].average.R, &stars[i].average.D, stars[i].measure.Xccd, stars[i].measure.Yccd, &image[0].coords);
+    stars[i].average.R = ohana_normalize_angle (stars[i].average.R);
   }
 

branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/get_image_patch.c

@@ -20 +20 @@
     XY_to_RD (&R, &D, Xpt[i]*image->NX, Ypt[i]*image->NY, &image->coords);
 
+    // fprintf (stderr, "%f %f -> %f %f\n", Xpt[i]*image->NX, Ypt[i]*image->NY, R, D); 
     Rmin = MIN(R,Rmin);
     Rmax = MAX(R,Rmax);
@@ -32 +33 @@
   region->Dmin = Dmin;
   region->Dmax = Dmax;
+
+  // fprintf (stderr, "%f %f , %f %f\n", Rmin, Rmax, Dmin, Dmax);
 
   return region;
@@ -48 +51 @@
 int SkyRegionsOverlap (SkyRegion *region, SkyRegion *patch) {
 
-  if (region->Rmin >= patch->Rmax) return FALSE;
-  if (region->Rmax <= patch->Rmin) return FALSE;
+  // what are the assumptions?  
+  // * patch->Rmin,Rmax are on the same side of the boundary 
+  //   * [-2,-1], [-2,+1], [358,361] : all OK
+  //   * [358, 2] : NOT OK
+
+  myAssert (patch->Rmin <= patch->Rmax, "impossible!");
+
   if (region->Dmin >= patch->Dmax) return FALSE;
   if (region->Dmax <= patch->Dmin) return FALSE;
+
+  int keepR0 = FALSE;
+  int keepR1 = FALSE;
+  int keepR2 = FALSE;
+  if (patch->Rmin < 0.0) {
+    float Rmin = patch->Rmin + 360.0;
+    float Rmax = patch->Rmax + 360.0;
+    keepR0 = (Rmin <= region->Rmax) && (Rmax >= region->Rmin);
+  }
+  if (patch->Rmax > 360.0) {
+    float Rmin = patch->Rmin - 360.0;
+    float Rmax = patch->Rmax - 360.0;
+    keepR1 = (Rmin <= region->Rmax) && (Rmax >= region->Rmin);
+  }
+  keepR2 = (patch->Rmin <= region->Rmax) && (patch->Rmax >= region->Rmin);
+
+  int keep = keepR0 || keepR1 || keepR2;
+
+  if (!keep) return FALSE;
   return TRUE;
 }

branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fake_images.c

@@ -71 +71 @@
       // gpc1-specific (assumes grizy map to 10101, etc)
       fakeImage[N].photcode = fakeImage[0].photcode + iy*10 + ix;
+
+      char *ptr = strstr (fakeImage[N].name, "[PHU]");
+      if (ptr) {
+        snprintf (ptr, 7, "[XY%d%d]", ix, iy);
+      }
       N++;
     }

branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fake_stars.c

@@ -5 +5 @@
   // patch is generous region around image, but limited ot this image
   SkyRegion *imagePatch = get_image_patch (image);
+
+  // fprintf (stderr, "image patch: %f %f , %f %f\n", imagePatch->Rmin, imagePatch->Rmax, imagePatch->Dmin, imagePatch->Dmax);
 
   int Nstars = 0;
@@ -14 +16 @@
     // skylist matches catalog
 
+    // fprintf (stderr, "try catalog: %f %f , %f %f : ", skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmin, skylist[0].regions[i][0].Dmax);
+
     // XXX check that catalog overlaps patch
-    if (!SkyRegionsOverlap(skylist[0].regions[i], imagePatch)) continue;
+    if (!SkyRegionsOverlap(skylist[0].regions[i], imagePatch)) { continue; }
+    // fprintf (stderr, "skip\n"); 
+    // fprintf (stderr, "keep\n");
 
     // generate fake measurements for this image

branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fake_stars_catalog.c

@@ -50 +50 @@
   StarPar *starpar = catalog->starpar;
   
+  double Rmid = 0.5*(imagePatch->Rmin + imagePatch->Rmax);
+
   int i;
   for (i = 0; i < catalog->Naverage; i++) {
 
     // true position from src catalog
-    if (!IN_PATCH(imagePatch, average[i].R, average[i].D)) continue;
+    double ra = ohana_normalize_angle_to_midpoint (average[i].R, Rmid);
+
+    // int show = !strcmp (image->name, "o6600g0127o.672655.cm.955522.smf[XY01]") && ((average[i].R < 1.0) || (average[i].R > 359.0));
+    int show = FALSE;
+
+    if (show) fprintf (stderr, "try: %f (%f) %f : %f %f , %f %f", average[i].R, ra, average[i].D, imagePatch->Rmin, imagePatch->Rmax, imagePatch->Dmin, imagePatch->Dmax);
+
+    if (!IN_PATCH(imagePatch, ra, average[i].D)) {
+      if (show) fprintf (stderr, "\n");
+      continue;
+    }
+    if (show) fprintf (stderr, " * \n");
 
     int nStar = average[i].starparOffset;
@@ -85 +98 @@
     // tzero, timeRef are in UNIX seconds, Toffset should be in years
     double Toffset = (image->tzero - timeRef) / 365.25 / 86400.0;
+    // Toffset = 0.0; // XXX TEST
 
     // uR,uD in linear (arcsec / yr)
@@ -100 +114 @@
     double Dobs = Dtru + dDoff;
 
+    Robs = ohana_normalize_angle (Robs);
+
     double X, Y;
     RD_to_XY (&X, &Y, Robs, Dobs, &image->coords);
+
+    if ((fabs(Robs) < 0.1) || (fabs(Robs - 360.0) < 0.1)) {
+      // fprintf (stderr, "%f %f : %f %f : %f %f\n", Rtru, Dtru, Robs, Dobs, X, Y);
+    }
+
     if (X < 0) continue;
     if (Y < 0) continue;
     if (X > image->NX) continue;
     if (Y > image->NY) continue;
+
+    if (!strcmp (image->name, "o6600g0127o.672655.cm.955522.smf[XY01]")) {
+      // fprintf (stderr, "keep: %f %f\n", average[i].R, average[i].D);
+    }
 
     stars[Nstars].measure.Xccd       = X;

branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/BoundaryTree.c

@@ -371 +371 @@
 int BoundaryTreeCellCoords (BoundaryTree *tree, int *zone, int *band, double ra, double dec) {
 
+  ra = ohana_normalize_angle (ra);
+
   // first, find the containing zone
 

branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/coordops.c

@@ -174 +174 @@
     *dec = DEG_RAD*delta;
 
-    *ra = ohana_normalize_angle (*ra);
+    // *ra = ohana_normalize_angle (*ra);
 
     return (TRUE);
@@ -207 +207 @@
     *dec = delta + coords[0].crval2;
 
-    *ra = ohana_normalize_angle (*ra);
+    // *ra = ohana_normalize_angle (*ra);
 
     return (TRUE);

branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/skyregion_ops.c

@@ -268 +268 @@
   }
 
-  if( Rmax > (Rmin+180.)) {
-    
-      double temp = Rmax;
-      Rmax = Rmin;
-      Rmin = temp;
-  }
-
   list = SkyListByBounds (table, depth, Rmin, Rmax, Dmin, Dmax);
 
@@ -283 +276 @@
 }
 
+// a region like -10 10 0 20 will become 350 10 0 20 and return all regions on both sides of 0,360
 SkyList *SkyListByBounds (SkyTable *table, int depth, double Rmin, double Rmax, double Dmin, double Dmax) {
 

branches/eam_branches/ipp-20140904/Ohana/src/opihi/dvo/objectcoverage.c

@@ -163 +163 @@
     for (xs = 0; xs < Nx; xs++) {
       status = XY_to_RD (&r, &d, (double)(xs), (double)(ys), &coords);
-      status &= (r >= 0);
-      status &= (r <= 360);
+      r = ohana_normalize_angle (r);
       if (status) {
         V[ys*Nx + xs] = ShowDensity ?  0 : 2;

branches/eam_branches/ipp-20140904/Ohana/src/opihi/dvo/skycoverage.c

@@ -227 +227 @@
     for (xs = 0; xs < Nx; xs++) {
       status = XY_to_RD (&r, &d, (double)(xs), (double)(ys), &coords);
-      status &= (r >= 0);
-      status &= (r <= 360);
+      r = ohana_normalize_angle (r);
+
       // are we in a part of the projection covering the sky or not?
       if (status) {

branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/UpdateObjects.c

@@ -357 +357 @@
     XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &coords);
     if (XVERB) fprintf (stderr, "project: %f %f : %f %f : %f\n", fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p);
+    if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
+    // XXX : does this make sense at 0,360 boundary?
 
     fitPM.p  = fitPM.dp  = 0.0;
@@ -408 +410 @@
     average[0].flags |= ID_STAR_FIT_PAR;
     fitStats->Npar ++;
+
+    if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
 
     // XXX a hard-wired hack...

branches/eam_branches/ipp-20140904/Ohana/src/relphot/src/MosaicOps.c

@@ -641 +641 @@
       Dmax = MAX (Dmax, D);
 
+      // XXX : this probably does not handle mosaics at RA = 0,360 well
+
       /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */