Index: trunk/magic/remove/src/diffedpixels.c
===================================================================
--- trunk/magic/remove/src/diffedpixels.c	(revision 27370)
+++ trunk/magic/remove/src/diffedpixels.c	(revision 27371)
@@ -8,4 +8,9 @@
 static int xLeft(psPlane *pt, int y);
 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
@@ -91,5 +96,5 @@
         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];
@@ -120,45 +125,58 @@
             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
-    fprintf(stderr, "\nSKYCELL %s Type: %s\n", fileName, type);
-    for (int i=0; i<4; i++) {
-        fprintf(stderr, "%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
-    // Now set the touched pixels
-    int ymin = fmax(0, pt[1].y );
-    int ymax = fmin(pt[3].y + 0.5, sf->diffedPixels->numRows - 1);
+
+    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
     for (int y = ymin ; y <= ymax; 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)
@@ -189,11 +207,16 @@
 {
     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);
-    }
+
+    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);
 }
@@ -202,90 +225,26 @@
 {
     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);
+    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);
 }
-
-/*
- * 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
-
-Type 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 >= pt2.y
-**************************
-
-                3
-              C
-        0----------------
-              B
-      ----------------2
-              A
-            1
-
-**************************
-If The left side corners have the same x we also have a Type 1
-
-        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
-
-*/
 
 static void
 nameCorners(psPlane pt[4])
 {
-    // sort points top to bottom
+    // sort points top to bottom (in case of ties take favor rightmost)
     int i;
     int j;
     for (i=0; i<4; i++) {
-        for (j = 0; j < 4; j++) {
-            if (i == j) continue;
-            if (pt[i].y > pt[j].y) {
+        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];
@@ -312,75 +271,7 @@
         tr = pt[0];
     }
-    int type = 0;
-    int type3 = 0;
-    if (round(tl.x) == round(bl.x)) {
-        type3 = 1;
-        if (tl.y <= tr.y) {
-            type = 1;
-        } else {
-            type = 2;
-        }
-    } else if (tl.y >= tr.y) {
-        type = 1;
-    } else {
-        type = 2;
-    }
-    if (type == 1) {
-        pt[0] = bl;
-        pt[1] = br;
-        pt[2] = tr;
-        pt[3] = tl;
-    } else {
-        pt[0] = tl;
-        pt[1] = bl;
-        pt[2] = br;
-        pt[3] = tr;
-    }
-
-    // Now check the results match our requirements
-    // pt3 is above pt0
-    if (pt[3].y <= pt[0].y) {
-        fprintf(stderr, "ERROR calculating diff overlap\n");
-        fprintf(stderr, "pt[3].y (%f) <= pt[0].y (%f)\n", pt[3].y, pt[0].y);
-        streaksExit("", PS_EXIT_PROG_ERROR);
-    }
-    // pt3 is above pt1
-    if (pt[3].y <= pt[1].y) {
-        fprintf(stderr, "ERROR calculating diff overlap\n");
-        fprintf(stderr, "pt[3].y (%f) <= pt[1].y (%f)\n", pt[3].y, pt[1].y);
-        streaksExit("", PS_EXIT_PROG_ERROR);
-    }
-    // pt2 is above pt1
-    if (pt[2].y <= pt[1].y) {
-        fprintf(stderr, "ERROR calculating diff overlap\n");
-        fprintf(stderr, "pt[2].y (%f) <= pt[1].y (%f)\n", pt[2].x, pt[1].x);
-        streaksExit("", PS_EXIT_PROG_ERROR);
-    }
-    // pt1 is to the right of pt0
-    if (pt[1].x <= pt[0].x) {
-        fprintf(stderr, "ERROR calculating diff overlap\n");
-        fprintf(stderr, "pt[1].x (%f) <= pt[0].x (%f)\n", pt[1].x, pt[0].x);
-        streaksExit("", PS_EXIT_PROG_ERROR);
-    }
-    // pt2 is to the right of pt0
-    if (pt[2].x <= pt[0].x) {
-        fprintf(stderr, "ERROR calculating diff overlap\n");
-        fprintf(stderr, "pt[2].x (%f) <= pt[0].x (%f)\n", pt[2].x, pt[0].x);
-        streaksExit("", PS_EXIT_PROG_ERROR);
-    }
-    // pt2 is to the right of pt3
-    if (pt[2].x <= pt[3].x) {
-        fprintf(stderr, "ERROR calculating diff overlap\n");
-        fprintf(stderr, "pt[2].x (%f) <= pt[3].x (%f)\n", pt[2].x, pt[3].x);
-        streaksExit("", PS_EXIT_PROG_ERROR);
-    }
-
-    // pt2 is below or at the same y as pt3
-    // This can happen with some strange shapes. Accept this with type3
-    // since our edge conditions are still satisfied
-    if (!type3 && (pt[2].y > pt[3].y)) {
-        fprintf(stderr, "ERROR calculating diff overlap\n");
-        fprintf(stderr, "pt[2].y (%f) > pt[3].y (%f)\n", pt[2].y, pt[3].y);
-        streaksExit("", PS_EXIT_PROG_ERROR);
-    }
-}
+    pt[BL] = bl;
+    pt[BR] = br;
+    pt[TL] = tl;
+    pt[TR] = tr;
+}