Changeset 37680

Timestamp:

Nov 27, 2014, 5:06:54 AM (12 years ago)

Author:

eugene

Message:

debugging grid fitting; adding pretty print

Location:

branches/eam_branches/ipp-20140904/Ohana/src/libdvo

Files:

• include/libdvo_astro.h (modified) (2 diffs)
• src/AstromOffsetMapIO.c (modified) (2 diffs)
• src/AstromOffsetMapOps.c (modified) (2 diffs)
• src/AstromOffsetMapUtils.c (modified) (6 diffs)

branches/eam_branches/ipp-20140904/Ohana/src/libdvo/include/libdvo_astro.h

@@ -64 +64 @@
   int Nmap;
   AstromOffsetMap **map;
-  int *IDtoSeq;
+   int *imageIDtoTableSeq;
   int MaxImageID;
   int MaxTableID;
@@ -174 +174 @@
 void AstromOffsetTableFree(AstromOffsetTable *table);
 int AstromOffsetTableAddMapFromImage (AstromOffsetTable *table, Image *image);
+void AstromOffsetMapPrint (AstromOffsetMap *map, char *filename);
+int AstromOffsetMapRepair (AstromOffsetMap *map, int xdir);
 
 # endif

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

@@ -116 +116 @@
 
   // generate the index and init values to -1
-  ALLOCATE (table->IDtoSeq, int, MaxImageID + 1);
+  ALLOCATE (table->imageIDtoTableSeq, int, MaxImageID + 1);
   for (i = 0; i <= MaxImageID; i++) {
-    table->IDtoSeq[i] = -1;
+    table->imageIDtoTableSeq[i] = -1;
   }
 
@@ -124 +124 @@
   for (i = 0; i < Nmap; i++) {
     int ImageID = map_disk[i].imageID;
-    myAssert (table->IDtoSeq[ImageID] == -1, "oops, duplicate image IDs");
-    table->IDtoSeq[ImageID] = i;
+    myAssert (table->imageIDtoTableSeq[ImageID] == -1, "oops, duplicate image IDs");
+    table->imageIDtoTableSeq[ImageID] = i;
   }
 

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

@@ -295 +295 @@
   ALLOCATE (Empty, int, Nx*Ny);
 
+  // find the max value
+  double MaxPivot = 0.0;
+  for (i = 0; i < Nx*Ny; i++) {
+    MaxPivot = MAX(A[i][i], MaxPivot);
+  }
+
+  // elements of A[i][j] are just the fractional weight of the number of input points to the pixel.
+  // any pixels which have < 5% contribution are not very well constrained.  NAN them for now...
+  double MinPivot = 0.05 * MaxPivot;
   for (i = 0; i < Nx*Ny; i++) {
     Empty[i] = 0;
-    if (fabs(A[i][i]) < 1.e-2) {
+    if (fabs(A[i][i]) < MinPivot) {
       Empty[i] = 1;
       for (j = 0; j < Nx*Ny; j++) {
@@ -352 +361 @@
 }
 
+// this function repairs an image with NAN pixels (only valid for a small-scale map -- no robust mean)
+int AstromOffsetMapRepair (AstromOffsetMap *map, int xdir) {
+
+  // we are going to repair the image by:
+  // 1) finding NAN pixels
+  // 2) if any of the neighbors are valid,
+  //    replace with the mean of the neighbors
+  // 3) otherwise, replace with the image mean
+
+  int Nx = map->Nx;
+  int Ny = map->Ny;
+
+  int ix, iy;
+
+  float **Vfix = NULL;
+  ALLOCATE (Vfix, float *, Nx);
+  for (ix = 0; ix < Nx; ix++) {
+    ALLOCATE (Vfix[ix], float, Ny);
+    for (iy = 0; iy < Ny; iy++) {
+      Vfix[ix][iy] = NAN;
+    }
+  }
+
+  // find the global mean
+  float mean = 0.0;
+  float npix = 0.0;
+
+  for (ix = 0; ix < Nx; ix++) {
+    for (iy = 0; iy < Ny; iy++) {
+      float value = xdir ? map->dXv[ix][iy] : map->dYv[ix][iy];
+      if (isnan(value)) continue;
+      mean += value;
+      npix += 1.0;
+    }
+  }
+  mean /= npix;
+
+  // find the NAN pixels:
+  for (ix = 0; ix < Nx; ix++) {
+    for (iy = 0; iy < Ny; iy++) {
+
+      float value = xdir ? map->dXv[ix][iy] : map->dYv[ix][iy];
+      if (!isnan(value)) {
+        Vfix[ix][iy] = value;
+        continue;
+      }
+
+      // find mean of all possible neighbors
+      float meanLocal = 0.0;
+      float npixLocal = 0.0;
+
+      int jx, jy;
+      for (jx = -1; jx <= +1; jx++) {
+        int nx = ix + jx;
+        if (nx < 0) continue;
+        if (nx >= Nx) continue;
+        for (jy = -1; jy <= +1; jy++) {
+          int ny = iy + jy;
+          if (ny < 0) continue;
+          if (ny >= Ny) continue;
+          float value = xdir ? map->dXv[nx][ny] : map->dYv[nx][ny];
+          if (isnan(value)) continue;
+
+          meanLocal += value;
+          npixLocal += 1.0;
+        }
+      }
+      // if there are no valid (non-NAN) local pixels, use global mean
+      meanLocal = (npixLocal > 0.0) ? meanLocal / npixLocal : mean;
+      Vfix[ix][iy] = meanLocal;
+    }
+  }
+    
+  // replace the bad pixels:
+  for (ix = 0; ix < Ny; ix++) {
+    for (iy = 0; iy < Ny; iy++) {
+      if (xdir) {
+        map->dXv[ix][iy] = Vfix[ix][iy];
+      } else {
+        map->dYv[ix][iy] = Vfix[ix][iy];
+      }
+    }
+  }
+  for (ix = 0; ix < Nx; ix++) {
+    free (Vfix[ix]);
+  }
+  free (Vfix);
+
+  return TRUE;
+}
+
 int dump_map_data (float *x, float *y, float *f, int Npts, char *filename) {
 

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

@@ -6 +6 @@
   // each of the maps so we can assign the map to the image.coords.offsetMap entry
 
-  // we have the lookup table of table->IDtoSeq to find the match by imageID
-  
-  // we are just going to assume that table->IDtoSeq is correct and complete
+  // we have the lookup table of table->imageIDtoTableSeq to find the match by imageID
+  
+  // we are just going to assume that table->imageIDtoTableSeq is correct and complete
   off_t i;
   for (i = 0; i < Nimages; i++) {
@@ -15 +15 @@
     if (imageID > table->MaxImageID) continue;
     
-    int seq = table->IDtoSeq[imageID];
+    int seq = table->imageIDtoTableSeq[imageID];
     if (seq < 0) continue;
     if (seq >= table->Nmap) continue; // this one is probably not valid, right?
@@ -83 +83 @@
     FREE (table->map[i][0].dYv);
   }
-  FREE (table->IDtoSeq);
+  FREE (table->imageIDtoTableSeq);
   FREE (table->map);
 }
@@ -101 +101 @@
     int oldMaxID = table->MaxImageID;
     table->MaxImageID = image->imageID;
-    REALLOCATE (table->IDtoSeq, int, table->MaxImageID + 1);
+    REALLOCATE (table->imageIDtoTableSeq, int, table->MaxImageID + 1);
     for (i = oldMaxID + 1; i < table->MaxImageID + 1; i++) {
-      table->IDtoSeq[i] = -1;
-    }
-  }
-  myAssert (table->IDtoSeq[image->imageID] == -1, "table IDtoSeq not initiazed or image collision");
-  table->IDtoSeq[image->imageID] = Nmap;
+      table->imageIDtoTableSeq[i] = -1;
+    }
+  }
+  int seq = table->imageIDtoTableSeq[image->imageID];
+  if (seq != -1) {
+    // already assigned (but not used)
+    // XXX this is probably wrong : I need to mark failed maps so they are not saved.
+    image[0].coords.offsetMap = table->map[seq];    
+    return TRUE;
+  }
+
+  myAssert (table->imageIDtoTableSeq[image->imageID] == -1, "table IDtoSeq not initiazed or image collision");
+  table->imageIDtoTableSeq[image->imageID] = Nmap;
   
   table->map[Nmap] = AstromOffsetMapInit (Nx, Ny);
@@ -129 +137 @@
   REALLOCATE (table->map, AstromOffsetMap *, table->Nmap);
 
-  // find the imageID and update the IDtoSeq allocation if needed
+  // find the imageID and update the imageIDtoTableSeq allocation if needed
   off_t i;
   if (image->imageID > table->MaxImageID) {
     int oldMaxID = table->MaxImageID;
     table->MaxImageID = image->imageID;
-    REALLOCATE (table->IDtoSeq, int, table->MaxImageID + 1);
+    REALLOCATE (table->imageIDtoTableSeq, int, table->MaxImageID + 1);
     for (i = oldMaxID + 1; i < table->MaxImageID + 1; i++) {
-      table->IDtoSeq[i] = -1;
-    }
-  }
-  myAssert (table->IDtoSeq[image->imageID] == -1, "table IDtoSeq not initiazed or image collision");
-  table->IDtoSeq[image->imageID] = Nmap;
+      table->imageIDtoTableSeq[i] = -1;
+    }
+  }
+  myAssert (table->imageIDtoTableSeq[image->imageID] == -1, "table IDtoSeq not initiazed or image collision");
+  table->imageIDtoTableSeq[image->imageID] = Nmap;
   
   table->map[Nmap] = image[0].coords.offsetMap;
@@ -158 +166 @@
 
   // generate the index and init values to -1
-  ALLOCATE (table->IDtoSeq, int, 1);
-  table->IDtoSeq[0] = -1;
+  ALLOCATE (table->imageIDtoTableSeq, int, 1);
+  table->imageIDtoTableSeq[0] = -1;
 
   return table;
 }
+
+void AstromOffsetMapPrint (AstromOffsetMap *map, char *filename) {
+
+  FILE *f = stderr;
+  if (strcasecmp(filename, "stderr") && strcasecmp(filename, "stdout")) {
+    f = fopen (filename, "w");
+    if (!f) {
+      fprintf (stderr, "failed to open output file %s\n", filename);
+    }
+  }
+
+  fprintf (f, "imageID: %d, tableID: %d (dX: %f, dY: %f)\n", map->imageID, map->tableID, map->dX, map->dY);
+  int ix, iy;
+  fprintf (f, "dXv map:\n");
+  for (ix = 0; ix < map->Nx; ix++) {
+    for (iy = 0; iy < map->Ny; iy++) {
+      fprintf (f, "%9.5f ", map->dXv[ix][iy]);
+    }
+    fprintf (f, "\n");
+  }
+  fprintf (f, "dYv map:\n");
+  for (ix = 0; ix < map->Nx; ix++) {
+    for (iy = 0; iy < map->Ny; iy++) {
+      fprintf (f, "%9.5f ", map->dYv[ix][iy]);
+    }
+    fprintf (f, "\n");
+  }
+  if (f != stderr) fclose (f);
+  return;
+}