Changeset 27517 for branches/eam_branches/20100225/magic/remove/src/diffedpixels.c

Timestamp:

Mar 29, 2010, 3:55:49 PM (16 years ago)

Author:

eugene

Message:

update merges from trunk

Location:

branches/eam_branches/20100225

Files:

• . (modified) (1 prop)
• magic/remove/src/diffedpixels.c (modified) (7 diffs)

branches/eam_branches/20100225/magic/remove/src/diffedpixels.c

@@ -2 +2 @@
 #include "pmAstrometryWCS.h"
 
-// #define DEBUG_PRINT 1
+#define DEBUG_PRINT 1
 
 static void addTouchedPixels(streakFiles *sf, psString fileName);
@@ -9 +9 @@
 static int xRight(psPlane *pt, int y);
 
+#define  BL 0
+#define  BR 1
+#define  TL 2
+#define  TR 3
+
 // Examine the set of diff skycells and compute the pixels that were
 // included in difference processing
@@ -37 +42 @@
         psString filename = psArrayGet(skycells, i);
         psString resolved_name = pmConfigConvertFilename(filename, config, false, false);
-
+        if (resolved_name == NULL) {
+            psError(PS_ERR_IO, false, "failed to resolve skycell file: %s", filename);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
         addTouchedPixels(sf, resolved_name);
         psFree(resolved_name);
@@ -88 +96 @@
         sf->diffedPixels->numCols, sf->diffedPixels->numRows);
 
-    // convert corners of skycell to sky coordinates
+    // convert corners of skycell to chip coordinates
     // compute overlap with the chip using astrometry
     psPlane pt[4];
@@ -117 +125 @@
             streaksExit("", PS_EXIT_DATA_ERROR);
         }
-        pt[i].x = p.x;
-        pt[i].y = p.y;
+        pt[i].x = round(p.x);
+        pt[i].y = round(p.y);
     }
 
     psFree(wcs);
 
-    // put the corners in the desired order (see comments below)
-
+    // Identify the bottom left, bottom right, top left, and top right corners
     nameCorners(pt);
-    psString type;
-    if (pt[0].y < pt[2].y ) {
-        type = "1";
-    } else {
-        type = "2";
-    }
+
 #ifdef DEBUG_PRINT
-    printf("\nSKYCELL %s Type: %s\n", fileName, type);
-    for (int i=0; i<4; i++) {
-        printf("%f %f\n", pt[i].x, pt[i].y);
-    }
+    fprintf(stderr, "\nSKYCELL %s\n", fileName);
+    fprintf(stderr, "BL %f %f\n", pt[BL].x, pt[BL].y);
+    fprintf(stderr, "BR %f %f\n", pt[BR].x, pt[BR].y);
+    fprintf(stderr, "TL %f %f\n", pt[TL].x, pt[TL].y);
+    fprintf(stderr, "TR %f %f\n", pt[TR].x, pt[TR].y);
 #endif
+
+    int xmax = sf->diffedPixels->numCols - 1;
+    int ymin;
+    int ymax; 
+    
+    if (pt[BL].y <= pt[BR].y) {
+        ymin = fmax(0, pt[BL].y);
+    } else {
+        ymin = fmax(0, pt[BR].y);
+    }
+    if (pt[TR].y >= pt[TL].y) {
+        ymax = fmin(pt[TR].y, sf->diffedPixels->numRows - 1);
+    } else {
+        ymax = fmin(pt[TL].y, sf->diffedPixels->numRows - 1);
+    }
+        
+#if (DEBUG_PRINT > 1)
+    fprintf(stderr, "\nxmin: %d xmax: %d\n", 0, xmax);
+    fprintf(stderr, "\nymin: %d ymax: %d\n", ymin, ymax);
+#endif
+
     // Now set the touched pixels
-    int ymin = fmax(0, pt[1].y );
-    int ymax = fmin(pt[3].y + 0.5, sf->diffedPixels->numRows - 1);
-#if (DEBUG_PRINT > 1)
-    printf("\nymin: %d ymax: %d\n", ymin, ymax);
-#endif
     for (int y = ymin ; y <= ymax; y++) {
-        int xleft  = xLeft(pt, y);
-        int xright = xRight(pt, y);
-
+
+        int xleft  = xLeft(pt, y) - 1;
+        int xright = xRight(pt, y) + 1;
+
+        if (xright < 0) {
+            continue;
+        }
+        if (xleft > xmax) {
+            continue;
+        }
         if (xleft < 0) {
             xleft = 0;
         }
-        if (xleft > sf->diffedPixels->numCols) {
-            continue;
-        }
-        if (xright < 0) {
-            continue;
-        }
-        if (xright >= sf->diffedPixels->numCols) {
-            xright = sf->diffedPixels->numCols - 1;
+        if (xright > xmax) {
+            xright = xmax;
         }
 #if (DEBUG_PRINT > 1)
-        printf("  y: %d xleft: %d xright: %d\n", y, xleft, xright);
+        fprintf(stderr, "  y: %d xleft: %d xright: %d\n", y, xleft, xright);
 #endif
 
@@ -187 +207 @@
 {
     double x_d;
-    if (y < pt[0].y) {
-        // left boundary is the line from pt 0 to pt 1
-        x_d = xOfY(&pt[0], &pt[1], y);
-    } else {
-        // left boundary is the line from pt 0 to pt 3
-        x_d = xOfY(&pt[0], &pt[3], y);
-    }
-    return (int) (x_d + 0.5);
+
+    if ((y <= pt[BL].y ) && ((pt[BL].y >= pt[BR].y))) {
+        // left boundary is the line from bottom left to bottom right
+        x_d = xOfY(&pt[BL], &pt[BR], y);
+    } else if ( y <= pt[TL].y ) {
+        // left boundary is the line from Bottom left to top left
+        x_d = xOfY(&pt[BL], &pt[TL], y);
+    } else {
+        // left boundary is the line from top left to top right
+        x_d= xOfY(&pt[TL], &pt[TR], y);
+    }
+
+    return (int) floor(x_d);
 }
 
@@ -200 +225 @@
 {
     double x_d;
-    if (y < pt[2].y) {
-        // right boundary is the line from pt 1 to pt 2
-        x_d = xOfY(&pt[1], &pt[2], y);
-
-    } else {
-        // right boundary is the line from pt 2 to pt 3
-        x_d = xOfY(&pt[2], &pt[3], y);
-    }
-    return (int) x_d;
-}
-
-/*
- * To compute the overlap of a quadrilateral with a set of
- * points we transform the 4 corners to image coordinates
- * and name the 4 points of the quad
-            pt 0 is left most (bottom most corner)
-            pt 1 is bottom most (right most corner)
-            pt 2 is right most (top most corner)
-            pt 3 is top most (left most corner)
-
-* Then we can divide the quad into 3 regions A B and C based
-* on the y coordinate
-
-* (in the degenerate case of a rectangle aligned to the axes
-* we only have 1 region)
-* In each of the three regions A, B, and C the test for whether a
-* point is contained in the quad on scanline y is simply
-
-        x >= xleft(y) x < xright(y)
-
-* where xleft(y) and xright(y) are the bounding lines at the given
-* y coordinate.
-
-The following diagrams show some examples. The 4 points are labeled 0 to 3
-The three regions are labeld A B and C, the boundaries between the regions
-ocur at pt0.y and pt2.y
-
-Example 1
-
-                3
-                 C
-         ---------------2           left boundary:  line 0_1 y < pt0.y
-                 B                                  line 0_3 y >= pt0.y
-          0--------------
-                 A                  right boundary: line 1_2 y < pt2.y
-                   1                                line 2_3  y >= pt.y
-**************************
-Example 2
-                3
-              C
-        0----------------
-              B
-      ----------------2
-              A
-            1
-
-**************************
-Example 3
-
-        3       2
-
-            B
-
-        0       1
-
-(region A and C are empty in the case where the points form a rectangle)
-*/
-
-/*
-    Name the corners
-    The following algorithm works for points that form a quadrilateral
-    I think it also works for the situation where 3 points are co-linear
-    and we have a triangle, but that isn't important for our purposes
-
-*/
+    if ((y <= pt[BR].y) && (pt[BR].y > pt[BL].y)) {
+        // right boundary is the line from Bottom right to Top right
+        x_d = xOfY(&pt[BL], &pt[BR], y);
+    } else if ( y<= pt[TR].y) {
+        // right boundary is the line from Top left to top right
+        x_d = xOfY(&pt[BR], &pt[TR], y);
+    } else {
+        x_d = xOfY(&pt[TL], &pt[TR], y);
+    }
+    return (int) ceil(x_d);
+}
 
 static void
 nameCorners(psPlane pt[4])
 {
-    psPlane   pt0, pt1, pt2, pt3;
-
-    /* find pt0 the left most (bottom most) corner */
-    int imin = 0;
-    int  min = (int) pt[0].x;
+    // sort points top to bottom (in case of ties take favor rightmost)
     int i;
-    for (i=1; i<4; i++) {
-        int x = pt[i].x;
-        if ((x < min) ||
-            ((x == min) && (pt[i].y < pt[imin].y))) {
-            imin = i;
-            min = x;
-        }
-    }
-    // remve pt0 from the list
-    pt0 = pt[imin];
-    for (i=imin; i<3; i++) {
-        pt[i] = pt[i+1];
-    }
-
-    // find the bottom most (right most) corner from the remaining 3
-    imin = 0;
-    min = pt[0].y;
-    for (i=1; i<3; i++) {
-        int y = pt[i].y;
-        if ((y < min) ||
-            ((y == min) && (pt[i].x > pt[imin].x)) ) {
-            imin = i;
-            min = y;
-        }
-    }
-    // remve pt1 from the list
-    pt1 = pt[imin];
-    for (i=imin; i<3; i++) {
-        pt[i] = pt[i+1];
-    }
-
-    // now find the right most (top most) of the remaining 2 points
-    if ((pt[0].x > pt[1].x) ||
-        ((pt[0].x == pt[1].x) && (pt[0].y > pt[1].y))) {
-        pt2 = pt[0];
-        pt3 = pt[1];
-    } else {
-        pt2 = pt[1];
-        pt3 = pt[0];
-    }
-    /* write the outputs */
-    pt[0] = pt0;
-    pt[1] = pt1;
-    pt[2] = pt2;
-    pt[3] = pt3;
-}
+    int j;
+    for (i=0; i<4; i++) {
+        for (j = i + 1; j < 4; j++) {
+            if ((pt[j].y > pt[i].y) || 
+                ((pt[j].y == pt[i].y) && (pt[j].x > pt[i].x))) {
+                psPlane tmpPt = pt[j];
+                pt[j] = pt[i];
+                pt[i] = tmpPt;
+            }
+        }
+    }
+    psPlane bl, br;
+    // sort bottom two points in x
+    if (pt[3].x < pt[2].x) {
+        bl = pt[3];
+        br = pt[2];
+    } else {
+        bl = pt[2];
+        br = pt[3];
+    }
+    // sort the top two points in x
+    psPlane tl, tr;
+    if (pt[0].x < pt[1].x) {
+        tl = pt[0];
+        tr = pt[1];
+    } else {
+        tl = pt[1];
+        tr = pt[0];
+    }
+    pt[BL] = bl;
+    pt[BR] = br;
+    pt[TL] = tl;
+    pt[TR] = tr;
+}