Changeset 14776

Timestamp:

Sep 7, 2007, 10:21:53 AM (19 years ago)

Author:

magnier

Message:

more work on making interpolation happy

Location:

branches/eam_branch_20070830/psLib/src/imageops

Files:

• psImagePixelInterpolate.c (modified) (11 diffs)
• psImagePixelInterpolate.h (modified) (3 diffs)

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

@@ -11 +11 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.1.2.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-03 20:29:07 $
+ *  @version $Revision: 1.1.2.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-07 20:21:51 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -38 +38 @@
 #include "psImageInterpolate.h"
 #include "psImagePixelInterpolate.h"
+
+#include "psFits.h"
+#include "psFitsImage.h"
 
 # define PS_IMAGE_ITER_STATE(XS,XE,YS,YE,N_MIN,TYPE) \
@@ -75 +78 @@
         for (int ix = 0; ix < mask->numCols; ix++) {
 
-            // skip the good pixels (unmasked)
-            if (!(mask->data.PS_TYPE_MASK_DATA[iy][ix] & maskVal)) { continue; }
+            // state of the good pixels (unmasked)
+            if (!(mask->data.PS_TYPE_MASK_DATA[iy][ix] & maskVal)) { 
+                // count good neighbor pixels (+ self)
+                int nGood = 0;
+                int minX = +1;
+                int maxX = -1;
+                int minY = +1;
+                int maxY = -1;
+                for (int jy = -1; jy <= +1; jy++) {
+                    /* stick to pixels in image grid */
+                    if (jy + iy < 0) { continue; }
+                    if (jy + iy >= mask->numRows) { continue; }
+                    for (int jx = -1; jx <= +1; jx++) {
+                        /* stick to pixels in image grid */
+                        if (jx + ix < 0) { continue; }
+                        if (jx + ix >= mask->numCols) { continue; }
+                        if (mask->data.PS_TYPE_MASK_DATA[iy+jy][ix+jx] & maskVal) { continue; }
+                        nGood ++;
+                        minX = PS_MIN (minX, jx);
+                        maxX = PS_MAX (maxX, jx);
+                        minY = PS_MIN (minY, jy);
+                        maxY = PS_MAX (maxY, jy);
+                    }
+                }
+                int dX = maxX - minX;
+                int dY = maxY - minY;
+                // what type of local interpolation can we use?
+                if ((nGood >= 6) && (dX == 2) && (dY == 2)) {
+                    result->data.S32[iy][ix] = PS_IMAGE_INTERPOLATE_GOOD2;
+                    continue; 
+                }
+                if (nGood >= 3) {
+                    result->data.S32[iy][ix] = PS_IMAGE_INTERPOLATE_GOOD1;
+                    continue; 
+                }
+                result->data.S32[iy][ix] = PS_IMAGE_INTERPOLATE_GOOD0;
+                continue; 
+            }
 
             // examine the neighbors.  If at least 6 of the 8 surrounding, or 3 of the 4 corner
@@ -85 +124 @@
             // check for poor pixels
             PS_IMAGE_ITER_STATE (-1,+1,-1,+1,6, PS_IMAGE_INTERPOLATE_CENTER);
-            PS_IMAGE_ITER_STATE (-1,+0,-1,+0,3, PS_IMAGE_INTERPOLATE_LL);
-            PS_IMAGE_ITER_STATE (-1,+0,+0,+1,3, PS_IMAGE_INTERPOLATE_UL);
-            PS_IMAGE_ITER_STATE (+0,+1,-1,+0,3, PS_IMAGE_INTERPOLATE_LR);
-            PS_IMAGE_ITER_STATE (+0,+1,+0,+1,3, PS_IMAGE_INTERPOLATE_UR);
+            PS_IMAGE_ITER_STATE (-1,+0,-1,+0,3, PS_IMAGE_INTERPOLATE_UR);
+            PS_IMAGE_ITER_STATE (-1,+0,+0,+1,3, PS_IMAGE_INTERPOLATE_LR);
+            PS_IMAGE_ITER_STATE (+0,+1,-1,+0,3, PS_IMAGE_INTERPOLATE_UL);
+            PS_IMAGE_ITER_STATE (+0,+1,+0,+1,3, PS_IMAGE_INTERPOLATE_LL);
 
             nBad ++;
@@ -122 +161 @@
 
             switch (state->data.S32[iy][ix]) {
-              case PS_IMAGE_INTERPOLATE_GOOD: 
+              case PS_IMAGE_INTERPOLATE_GOOD0: 
+              case PS_IMAGE_INTERPOLATE_GOOD1: 
+              case PS_IMAGE_INTERPOLATE_GOOD2: 
                 // skip the good pixels
                 break; 
@@ -167 +208 @@
               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)
-                  image->data.F32[iy][ix] = image->data.F32[iy+1][ix+0] - image->data.F32[iy+0][ix+1] - image->data.F32[iy+1][ix+1];
+                  image->data.F32[iy][ix] = image->data.F32[iy+1][ix+1] - image->data.F32[iy+0][ix+1] - image->data.F32[iy+1][ix+0];
                   break; }
 
               case PS_IMAGE_INTERPOLATE_LR: {
                   // fit a plane to the 3 pixels at (0,1),(-1,0),(-1,1), extend to pixel at (0,0)
-                  image->data.F32[iy][ix] = image->data.F32[iy+1][ix+0] - image->data.F32[iy+0][ix-1] - image->data.F32[iy+1][ix-1];
+                  image->data.F32[iy][ix] = image->data.F32[iy+1][ix-1] - image->data.F32[iy+0][ix-1] - image->data.F32[iy+1][ix+0];
                   break; }
 
               case PS_IMAGE_INTERPOLATE_UL: {
                   // fit a plane to the 3 pixels at (0,-1),(1,0),(1,-1), extend to pixel at (0,0)
-                  image->data.F32[iy][ix] = image->data.F32[iy-1][ix+0] - image->data.F32[iy+0][ix+1] - image->data.F32[iy-1][ix+1];
+                  image->data.F32[iy][ix] = image->data.F32[iy-1][ix+1] - image->data.F32[iy+0][ix+1] - image->data.F32[iy-1][ix+0];
                   break; }
 
               case PS_IMAGE_INTERPOLATE_UR: {
                   // fit a plane to the 3 pixels at (0,-1),(-1,0),(-1,-1), extend to pixel at (0,0)
-                  image->data.F32[iy][ix] = image->data.F32[iy-1][ix+0] - image->data.F32[iy+0][ix-1] - image->data.F32[iy-1][ix-1];
+                  image->data.F32[iy][ix] = image->data.F32[iy-1][ix-1] - image->data.F32[iy+0][ix-1] - image->data.F32[iy-1][ix+0];
                   break; }
 
@@ -207 +248 @@
     assert (value->numRows == mask->numRows);
 
+# if (0)
+    psFits *fits = NULL;
+
+    fits = psFitsOpen ("xcoords.fits", "w");
+    psFitsWriteImage (fits, NULL, xCoord, 0, NULL);
+    psFitsClose (fits);
+
+    fits = psFitsOpen ("ycoords.fits", "w");
+    psFitsWriteImage (fits, NULL, yCoord, 0, NULL);
+    psFitsClose (fits);
+
+    fits = psFitsOpen ("value.fits", "w");
+    psFitsWriteImage (fits, NULL, value, 0, NULL);
+    psFitsClose (fits);
+# endif
+
     psImage *output = psImageAlloc (value->numCols, value->numRows, PS_TYPE_F32);
 
@@ -217 +274 @@
 
     // 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;
+    psPolynomial2D *poly2o = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2);
+    poly2o->mask[2][2] = 1;
+    poly2o->mask[2][1] = 1;
+    poly2o->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;
+    psPolynomial2D *poly1o = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 1, 1);
+    poly2o->mask[1][1] = 1;
 
     for (int iy = 0; iy < state->numRows; iy++) {
@@ -256 +313 @@
                   // df->n = n;
                   // psVectorFitPolynomial2D (poly, NULL, 0xff, f, df, x, y);
-                  psVectorFitPolynomial2D (poly2D, NULL, 0xff, f, NULL, x, y);
+                  psVectorFitPolynomial2D (poly2o, 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);
+                  // center of pixel is 0.5,0.5
+                  output->data.F32[iy][ix] = psPolynomial2DEval (poly2o, ix + 0.5, iy + 0.5);
                   break; }
 
@@ -286 +344 @@
                   f->n = n;
                   // df->n = n;
-                  // psVectorFitPolynomial2D (poly1D, NULL, 0xff, f, df, x, y);
-                  psVectorFitPolynomial2D (poly1D, NULL, 0xff, f, NULL, x, y);
+                  // psVectorFitPolynomial2D (poly1o, NULL, 0xff, f, df, x, y);
+                  psVectorFitPolynomial2D (poly1o, 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);
+                  output->data.F32[iy][ix] = psPolynomial2DEval (poly1o, ix + 0.5, iy + 0.5);
                   break; }
 
               case PS_IMAGE_INTERPOLATE_GOOD0: {
-                  output->data.F32[iy][ix] = value->data.F32[iy+jy][ix+jx];
+                  output->data.F32[iy][ix] = value->data.F32[iy][ix];
                   break; }
 
               default:
-                psAbort("impossible case in __func__");
+                // skip poor or bad pixels (interpolate later)
+                break;
             }
         }           
@@ -313 +372 @@
     psFree (f);
 
-    psFree (poly2D);
-    psFree (poly1D);
+    psFree (poly2o);
+    psFree (poly1o);
 
     psFree (output);

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

@@ -7 +7 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-02 20:29:50 $
+ *  @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-07 20:21:53 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -21 +21 @@
 // XXX make these all bit values?
 typedef enum {
-    PS_IMAGE_INTERPOLATE_GOOD   = 0x00,
+    PS_IMAGE_INTERPOLATE_GOOD   = 0x10,
+    PS_IMAGE_INTERPOLATE_GOOD0  = 0x11,
+    PS_IMAGE_INTERPOLATE_GOOD1  = 0x12,
+    PS_IMAGE_INTERPOLATE_GOOD2  = 0x13,
     PS_IMAGE_INTERPOLATE_BAD    = 0x01,
     PS_IMAGE_INTERPOLATE_CENTER = 0x02,
@@ -33 +36 @@
 psImage *psImagePixelInterpolateState (int *nBad, int *nPoor, psImage *mask, psMaskType maskVal);
 bool psImagePixelInterpolatePoor (psImage *image, psImage *state, psImage *mask, psMaskType maskVal);
+bool psImagePixelInterpolateCenter (psImage *value, psImage *xCoord, psImage *yCoord, psImage *state, psImage *mask, psMaskType maskVal);
 
 /// @}