Changeset 27371 for trunk/magic/remove/src/diffedpixels.c

Timestamp:

Mar 19, 2010, 6:32:46 PM (16 years ago)

Author:

bills

Message:

A skycell boundary transformed to chip coordinates can have 3 left edges. Deal with this by
changing the way we identify the corners to work with upper left, upper right, lower left and lower right
corners. When calculating the left and right boundaries deal up to 2 transitions

File:

• trunk/magic/remove/src/diffedpixels.c (modified) (6 diffs)

trunk/magic/remove/src/diffedpixels.c

@@ -8 +8 @@
 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 +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];
@@ -120 +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
-    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 +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);
-    }
+
+    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 +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);
+    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 +271 @@
         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;
+}