Mercurial > ldcheck / changeset

--- a/geometry.py	Sun Jan 21 15:03:38 2018 +0200
+++ b/geometry.py	Mon Jan 22 12:22:01 2018 +0200
@@ -318,6 +318,16 @@
             sum(self[i, k] * other[k, j] for k in range(3))
             for i, j in matrix_indices()
         ])
+    def __add__(self, other):
+        return Matrix3x3([
+            a + b
+            for a, b in zip(self.values, other.values)
+        ])
+    def __mul__(self, scalar):
+        return Matrix3x3([
+            x * scalar
+            for x in self.values
+        ])
 
 def complete_matrix(matrix, anchor):
     '''
@@ -349,3 +359,58 @@
             + transformation_matrix[10] * vertex.z \
             + transformation_matrix[11],
     )
+
+def transform_to_xy(polygon):
+    '''
+        Transforms the provided polygon into the XY plane. The polygon is
+        assumed to be planar.
+
+        Implementation is based on this StackOverflow answer:
+        https://math.stackexchange.com/a/476311
+    '''
+    v1, v2, v3 = polygon.vertices[:3]
+    a, b = v3 - v2, v1 - v2
+    normal = cross_product(a, b).normalized()
+    v = cross_product(normal, Vector(0, 0, 1))
+    cosine = dot_product(normal, Vector(0, 0, 1))
+    v_matrix = Matrix3x3([
+        0, -v.z, v.y,
+        v.z, 0, -v.x,
+        -v.y, v.x, 0,
+    ])
+    try:
+        rotation_matrix = Matrix3x3()
+        rotation_matrix += v_matrix
+        rotation_matrix += (v_matrix @ v_matrix) * (1 / (1 + cosine))
+    except ZeroDivisionError:
+        rotation_matrix = Matrix3x3()
+    full_matrix = complete_matrix(rotation_matrix, Vertex(0, 0, 0))
+    return Polygon([
+        transform(vertex = vertex, transformation_matrix = full_matrix)
+        for vertex in polygon.vertices
+    ])
+
+def vector_angle(vec_1, vec_2, normalized = False):
+    from math import acos, pi as π
+    cosine = dot_product(vec_1, vec_2)
+    try:
+        cosine /= vec_1.length() * vec_2.length()
+    except ZeroDivisionError:
+        return 0
+    angle = acos(cosine)
+    if normalized and angle > π / 2:
+        angle = π - angle
+    return angle
+
+def split_quadrilateral(polygon):
+    assert(len(polygon.vertices) == 4)
+    vertices = IndexRing(polygon.vertices)
+    for i in (0, 1):
+        a, b = vertices[0 + i], vertices[1 + i]
+        c, d = vertices[2 + i], vertices[3 + i]
+        yield Polygon([a, b, c]), Polygon([b, c, d])
+
+def triangle_plane_normal(polygon):
+    assert(len(polygon.vertices) == 3)
+    a, b, c = polygon.vertices[:3]
+    return cross_product(b - a, b - c)

--- a/tests/concave.py	Sun Jan 21 15:03:38 2018 +0200
+++ b/tests/concave.py	Mon Jan 22 12:22:01 2018 +0200
@@ -6,57 +6,31 @@
     # Returns whether all elements in container have the same sign
     return min(container) * max(container) >= 0
 
-def transform_to_xy(geometry):
-    a, b, c = geometry.vertices[:3]
-    
-
 def concave_test(model):
     for quadrilateral in model.quadrilaterals:
-        print([cross_product(v2 - v1, v3 - v1) for v1, v2, v3 in pairs(quadrilateral.geometry.vertices, count = 3)])
+        geometry = transform_to_xy(quadrilateral.geometry)
         z_scores = [
             cross_product(v2 - v1, v3 - v1).z
-            for v1, v2, v3 in pairs(quadrilateral.geometry.vertices, count = 3)
+            for v1, v2, v3 in pairs(geometry.vertices, count = 3)
         ]
-        print(z_scores)
         if not sign_consistency(z_scores):
-            yield warning(quadrilateral, 'Concave quadrilateral')
+            yield error(quadrilateral, 'Concave quadrilateral')
 
 def bowtie_quadrilateral_test(model):
-    for quadrilateral in model.quadrilaterals:
-        vertices = IndexRing(quadrilateral.geometry.vertices)
-        for i in (0, 1):
-            line1 = LineSegment(vertices[0 + i], vertices[1 + i])
-            line2 = LineSegment(vertices[2 + i], vertices[3 + i])
-            try:
-                line_intersection(line1, line2)
-            except NoIntersection:
-                pass
-            else:
-                yield error(quadrilateral, 'Bowtie quadrilateral')
-                break
-
-def vector_angle(vec_1, vec_2, normalized = False):
-    cosine = dot_product(vec_1, vec_2)
-    try:
-        cosine /= vec_1.length() * vec_2.length()
-    except ZeroDivisionError:
-        return 0
-    angle = acos(cosine)
-    if normalized and angle > π / 2:
-        angle = π - angle
-    return angle
+    ...
 
 def skew_test(model):
     for quadrilateral in model.quadrilaterals:
-        vertices = IndexRing(quadrilateral.geometry.vertices)
-        for i in (0, 1):
-            a, b = vertices[0 + i], vertices[1 + i]
-            c, d = vertices[2 + i], vertices[3 + i]
-            plane_1 = cross_product(b - a, d - a)
-            plane_2 = cross_product(d - c, b - c)
+        for triangles in split_quadrilateral(quadrilateral.geometry):
+            plane_1 = triangle_plane_normal(triangles[0])
+            plane_2 = triangle_plane_normal(triangles[1])
             angle = vector_angle(plane_1, plane_2, normalized = True)
             if angle > radians(0.1):
                 yield error(quadrilateral,
                     'Skew quadrilateral (plane angle {}°)',
                     '%.3f' % degrees(angle))
                 break
+
+manifest = {
+    'tests': {skew_test, concave_test},
+}