[Plplot-cvs] SF.net SVN: plplot:[10728] trunk/src/plfill.c

[<ai...@us...>]

Revision: 10728
          http://plplot.svn.sourceforge.net/plplot/?rev=10728&view=rev
Author:   airwin
Date:     2009-12-22 19:15:21 +0000 (Tue, 22 Dec 2009)

Log Message:
-----------
Upgrade positive_orientation logic from a simple rule which only works
consistently for convex polygons to a general rule which works for all
non-intersecting polygons whether convex or not.  By chance, this
-DUSE_FILL_INTERSECTION_POLYGON=ON bug fix makes no difference in the
results for example 25, but it will make a difference in the general case.

Modified Paths:
--------------
    trunk/src/plfill.c

Modified: trunk/src/plfill.c
===================================================================
--- trunk/src/plfill.c	2009-12-22 02:31:06 UTC (rev 10727)
+++ trunk/src/plfill.c	2009-12-22 19:15:21 UTC (rev 10728)
@@ -1930,30 +1930,40 @@
 }
 
 /* Decide if polygon has a positive orientation or not.
- * See http://en.wikipedia.org/wiki/Curve_orientation for details
- * of this simple determinate method.  */
+ * Note the simple algorithm given in
+ * http://en.wikipedia.org/wiki/Curve_orientation is incorrect for
+ * non-convex polygons.  Instead, for the general nonintersecting case
+ * use the polygon area method given by
+ * http://local.wasp.uwa.edu.au/~pbourke/geometry/polyarea/ where
+ * you project each edge of the polygon down to the X axis and take the
+ * area of the enclosed trapezoid.  The trapezoid areas outside the
+ * polygon are completely cancelled if you sum over all edges.  Furthermore,
+ * the sum of the trapezoid areas have terms which are zero when calculated
+ * with the telescoping rule so the final formula is quite simple. */
 int
 positive_orientation( PLINT n, const PLINT *x, const PLINT *y )
 {
-    PLFLT xa, ya, xb, yb, xc, yc, det;
+    PLINT i, im1;
+    /* Use PLFLT for all calculations to avoid integer overflows.  Also,
+     * the normal PLFLT has 64-bits which means you get exact integer
+     * arithmetic well beyond the normal integer overflow limits. */
+    PLFLT twice_area = 0.;
     if ( n < 3 )
     {
         plwarn( "positive_orientation: internal logic error, n < 3" );
         return 0;
     }
-    /* Use floating point to avoid integer overflows. */
-    xa  = x[0];
-    xb  = x[1];
-    xc  = x[2];
-    ya  = y[0];
-    yb  = y[1];
-    yc  = y[2];
-    det = ( xa * yb + xb * yc + xc * ya ) - ( xa * yc + xb * ya + xc * yb );
-    if ( det == 0. )
+    im1 = n - 1;
+    for ( i = 0; i < n; i++ )
     {
-        plwarn( "positive_orientation: internal logic error, det == 0." );
+        twice_area += (PLFLT) x[im1] * (PLFLT) y[i] - (PLFLT) x[i] * (PLFLT) y[im1];
+        im1         = i;
+    }
+    if ( twice_area == 0. )
+    {
+        plwarn( "positive_orientation: internal logic error, twice_area == 0." );
         return 0;
     }
     else
-        return det > 0.;
+        return twice_area > 0.;
 }


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.