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

Timestamp:

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

Author:

eugene

Message:

debugging grid fitting; adding pretty print

File:

• branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/AstromOffsetMapOps.c (modified) (2 diffs)

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) {