Fix methods Domain._sort_polygon_2d(), Domain._sort_polygon_3d()

diff --git a/pypol/polyhedra.py b/pypol/polyhedra.py
index 37f16e0..69ed2b2 100644 (file)
--- a/pypol/polyhedra.py
+++ b/pypol/polyhedra.py
@@ -5,7 +5,7 @@ import numbers
 from . import islhelper
 
 from .islhelper import mainctx, libisl
-from .geometry import GeometricObject, Point
+from .geometry import GeometricObject, Point, Vector
 from .linexprs import Expression, Symbol, Rational
 from .domains import Domain
 
@@ -206,16 +206,26 @@ class Polyhedron(Domain):
             constraints.append(sympy.Ge(inequality.tosympy(), 0))
         return sympy.And(*constraints)
 
+    @classmethod
+    def _polygon_inner_point(cls, points):
+        symbols = points[0].symbols
+        coordinates = {symbol: 0 for symbol in symbols}
+        for point in points:
+            for symbol, coordinate in point.coordinates():
+                coordinates[symbol] += coordinate
+        for symbol in symbols:
+            coordinates[symbol] /= len(points)
+        return Point(coordinates)
+
     @classmethod
     def _sort_polygon_2d(cls, points):
         if len(points) <= 3:
             return points
-        o = sum((Vector(point) for point in points)) / len(points)
-        o = Point(o.coordinates())
+        o = cls._polygon_inner_point(points)
         angles = {}
         for m in points:
             om = Vector(o, m)
-            dx, dy = (coordinate for symbol, coordinates in om.coordinates())
+            dx, dy = (coordinate for symbol, coordinate in om.coordinates())
             angle = math.atan2(dy, dx)
             angles[m] = angle
         return sorted(points, key=angles.get)
@@ -224,13 +234,18 @@ class Polyhedron(Domain):
     def _sort_polygon_3d(cls, points):
         if len(points) <= 3:
             return points
-        o = sum((Vector(point) for point in points)) / len(points)
-        o = Point(o.coordinates())
-        a, b = points[:2]
+        o = cls._polygon_inner_point(points)
+        a = points[0]
         oa = Vector(o, a)
-        ob = Vector(o, b)
         norm_oa = oa.norm()
-        u = (oa.cross(ob)).asunit()
+        for b in points[1:]:
+            ob = Vector(o, b)
+            u = oa.cross(ob)
+            if not u.isnull():
+                u = u.asunit()
+                break
+        else:
+            raise ValueError('degenerate polygon')
         angles = {a: 0.}
         for m in points[1:]:
             om = Vector(o, m)