Index: trunk/magic/remove/src/diffedpixels.c
===================================================================
--- trunk/magic/remove/src/diffedpixels.c	(revision 27102)
+++ trunk/magic/remove/src/diffedpixels.c	(revision 27223)
@@ -2,5 +2,5 @@
 #include "pmAstrometryWCS.h"
 
-// #define DEBUG_PRINT 1
+#define DEBUG_PRINT 1
 
 static void addTouchedPixels(streakFiles *sf, psString fileName);
@@ -136,7 +136,7 @@
     }
 #ifdef DEBUG_PRINT
-    printf("\nSKYCELL %s Type: %s\n", fileName, type);
+    fprintf(stderr, "\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, "%f %f\n", pt[i].x, pt[i].y);
     }
 #endif
@@ -145,5 +145,5 @@
     int ymax = fmin(pt[3].y + 0.5, sf->diffedPixels->numRows - 1);
 #if (DEBUG_PRINT > 1)
-    printf("\nymin: %d ymax: %d\n", ymin, ymax);
+    fprintf(stderr, "\nymin: %d ymax: %d\n", ymin, ymax);
 #endif
     for (int y = ymin ; y <= ymax; y++) {
@@ -164,5 +164,5 @@
         }
 #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
 
@@ -284,52 +284,93 @@
     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
     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];
+    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) {
+                psPlane tmpPt = pt[j];
+                pt[j] = pt[i];
+                pt[i] = tmpPt;
+            }
+        }
+    }
+    // sort bottom two points in x
+    if (pt[3].x < pt[2].x) {
+        psPlane tmpPt = pt[3];
+        pt[3] = pt[2];
+        pt[2] = tmpPt;
+    }
+    // now we know that p[3] point is pt1 - the lower right corner
+    pt1 = pt[3];
+    
+    // And pt[2] is left most of the bottom 2  (either pt0 or pt2)
+
+    // find which of the top two is to the left
+    if (pt[0].x < pt[1].x) {
+        pt3 = pt[0];
+        // now decide which of pt[1] and pt[2] is pt0 - on the left)
+        if (pt[1].x < pt[2].x) {
+            pt0 = pt[1];
+            pt2 = pt[2];
+        } else {
+            pt0 = pt[2];
+            pt2 = pt[1];
+        }
+    } else {
         pt3 = pt[1];
-    } else {
-        pt2 = pt[1];
-        pt3 = pt[0];
-    }
-    /* write the outputs */
+        // now decide which of pt[0] and pt[2] is pt0 - on the left)
+        if (pt[0].x < pt[2].x) {
+            pt0 = pt[0];
+            pt2 = pt[2];
+        } else {
+            pt0 = pt[2];
+            pt2 = pt[0];
+        }
+    }
+
+    /* now write the outputs back to the input array*/
     pt[0] = pt0;
     pt[1] = pt1;
     pt[2] = pt2;
     pt[3] = pt3;
-}
+
+    // And Check the results
+    // y3 > y0
+    // y2 > y0
+    // y3 > y1
+    // x1 > x0
+    // x2 > x0
+    // x2 > x3
+    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);
+    }
+    if (pt[2].y <= pt[0].y) {
+        fprintf(stderr, "ERROR calculating diff overlap\n");
+        fprintf(stderr, "pt[2].y (%f) <= pt[0].y (%f)\n", pt[2].y, pt[0].y);
+        streaksExit("", PS_EXIT_PROG_ERROR);
+    }
+    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);
+    }
+    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);
+    }
+    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);
+    }
+    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);
+    }
+}