Changeset 14726 for branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.c

Timestamp:

Sep 3, 2007, 10:29:07 AM (19 years ago)

Author:

magnier

Message:

updates to image map code, tests

File:

• branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.c (modified) (3 diffs)

branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.c

@@ -11 +11 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.1.2.4 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-02 20:29:50 $
+ *  @version $Revision: 1.1.2.5 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-03 20:29:07 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -163 +163 @@
                   break; }
 
+                // XXX should I use 1 1D polynomial fitting all unmasked pixels in the 3x3 grid?
+                // XXX that would automatically extend to regions where only 2 pixels are valid...
               case PS_IMAGE_INTERPOLATE_LL: {
                   // fit a plane to the 3 pixels at (0,1),(1,0),(1,1), extend to pixel at (0,0)
@@ -196 +198 @@
     return true;
 }
-
     
+// interpolate the good pixels to their true centers
+bool psImagePixelInterpolateCenter (psImage *value, psImage *xCoord, psImage *yCoord, psImage *state, psImage *mask, psMaskType maskVal) {
+
+    assert (value->numCols == state->numCols);
+    assert (value->numRows == state->numRows);
+    assert (value->numCols == mask->numCols);
+    assert (value->numRows == mask->numRows);
+
+    psImage *output = psImageAlloc (value->numCols, value->numRows, PS_TYPE_F32);
+
+    // allocate the vectors for the 2nd order fit below
+    psVector *f  = psVectorAlloc (9, PS_TYPE_F32);
+    // XXX if we add the weight above, include df
+    // psVector *df = psVectorAlloc (9, PS_TYPE_F32); 
+    psVector *x  = psVectorAlloc (9, PS_TYPE_F32);
+    psVector *y  = psVectorAlloc (9, PS_TYPE_F32);
+
+    // allocate a 2D polynomial to fit a quadratic to the valid neighbor pixels.
+    psPolynomial2D *poly2D = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2);
+    poly2D->mask[2][2] = 1;
+    poly2D->mask[2][1] = 1;
+    poly2D->mask[1][2] = 1;
+
+    // allocate a 2D polynomial to fit a plane to the valid neighbor pixels.
+    psPolynomial2D *poly1D = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 1, 1);
+    poly2D->mask[1][1] = 1;
+
+    for (int iy = 0; iy < state->numRows; iy++) {
+        for (int ix = 0; ix < state->numCols; ix++) {
+
+            switch (state->data.S32[iy][ix]) {
+              case PS_IMAGE_INTERPOLATE_GOOD2: {
+                  // XXX is there a fit-image-region function?
+                  int n = 0;
+                  for (int jy = -1; jy <= +1; jy++) {
+                      // skip invalid pixels 
+                      if (jy + iy < 0) { continue; }
+                      if (jy + iy >= value->numRows) { continue; }
+                      for (int jx = -1; jx <= +1; jx++) {
+                          // skip invalid pixels 
+                          if (jx + ix < 0) { continue; } 
+                          if (jx + ix >= value->numCols) { continue; } 
+                          // skip masked pixels
+                          if (mask->data.PS_TYPE_MASK_DATA[iy+jy][ix+jx] & maskVal) { continue; }
+                          x->data.F32[n] = xCoord->data.F32[iy+jy][ix+jx];
+                          y->data.F32[n] = yCoord->data.F32[iy+jy][ix+jx];
+                          f->data.F32[n] = value->data.F32[iy+jy][ix+jx];
+                          // df->data.F32[n] = weight->data.F32[iy+jy][ix+jx];
+                          n++;
+                      }
+                  }
+                  // set vector lengths here
+                  x->n = n;
+                  y->n = n;
+                  f->n = n;
+                  // df->n = n;
+                  // psVectorFitPolynomial2D (poly, NULL, 0xff, f, df, x, y);
+                  psVectorFitPolynomial2D (poly2D, NULL, 0xff, f, NULL, x, y);
+                  // apply the fitted quadratic to get the poor pixel value
+                  output->data.F32[iy][ix] = psPolynomial2DEval (poly2D, ix, iy);
+                  break; }
+
+              case PS_IMAGE_INTERPOLATE_GOOD1: {
+                  // XXX is there a fit-image-region function?
+                  int n = 0;
+                  for (int jy = -1; jy <= +1; jy++) {
+                      // skip invalid pixels 
+                      if (jy + iy < 0) { continue; }
+                      if (jy + iy >= value->numRows) { continue; }
+                      for (int jx = -1; jx <= +1; jx++) {
+                          // skip invalid pixels 
+                          if (jx + ix < 0) { continue; } 
+                          if (jx + ix >= value->numCols) { continue; } 
+                          // skip masked pixels
+                          if (mask->data.PS_TYPE_MASK_DATA[iy+jy][ix+jx] & maskVal) { continue; }
+                          x->data.F32[n] = xCoord->data.F32[iy+jy][ix+jx];
+                          y->data.F32[n] = yCoord->data.F32[iy+jy][ix+jx];
+                          f->data.F32[n] = value->data.F32[iy+jy][ix+jx];
+                          // df->data.F32[n] = weight->data.F32[iy+jy][ix+jx];
+                          n++;
+                      }
+                  }
+                  // set vector lengths here
+                  x->n = n;
+                  y->n = n;
+                  f->n = n;
+                  // df->n = n;
+                  // psVectorFitPolynomial2D (poly1D, NULL, 0xff, f, df, x, y);
+                  psVectorFitPolynomial2D (poly1D, NULL, 0xff, f, NULL, x, y);
+                  // apply the fitted quadratic to get the poor pixel value
+                  output->data.F32[iy][ix] = psPolynomial2DEval (poly1D, ix, iy);
+                  break; }
+
+              case PS_IMAGE_INTERPOLATE_GOOD0: {
+                  output->data.F32[iy][ix] = value->data.F32[iy+jy][ix+jx];
+                  break; }
+
+              default:
+                psAbort("impossible case in __func__");
+            }
+        }           
+    }
+
+    for (int iy = 0; iy < value->numRows; iy++) {
+        for (int ix = 0; ix < value->numCols; ix++) {
+          if (mask->data.PS_TYPE_MASK_DATA[iy][ix] & maskVal) { continue; }
+          value->data.F32[iy][ix] = output->data.F32[iy][ix];
+        }
+    }
+
+    psFree (x);
+    psFree (y);
+    psFree (f);
+
+    psFree (poly2D);
+    psFree (poly1D);
+
+    psFree (output);
+    return true;
+}