add function ezdxf.math.is_vertex_order_ccw_3d()

docs/source/math/core.rst

@@ -185,6 +185,7 @@ Example for a closed collinear shape, which creates 2 additional vertices and th
     intersection_ray_polygon_3d
     intersection_ray_ray_3d
     is_planar_face
+    is_vertex_order_ccw_3d
     linear_vertex_spacing
     local_cubic_bspline_interpolation
     normal_vector_3p
@@ -251,6 +252,8 @@ Example for a closed collinear shape, which creates 2 additional vertices and th
 
 .. autofunction:: is_planar_face
 
+.. autofunction:: is_vertex_order_ccw_3d
+
 .. autofunction:: linear_vertex_spacing
 
 .. autofunction:: local_cubic_bspline_interpolation

notes/pages/CHANGELOG.md

@@ -9,6 +9,7 @@
 	- NEW: `BSpline.measure()` method
 	- NEW: `BSpline.split()` method
 	- NEW: `BSpline.insert_knot()` and `BSpline.knot_refinement()` supports rational splines
+	- NEW: function `ezdxf.math.is_vertex_order_ccw_3d()`
 	- BUGFIX: Exported `MESH` entities  without vertices or faces create invalid DXF files
 		- {{issue 1219}}
 	- BUGFIX: `pickle` support added by #mbway

src/ezdxf/math/construct3d.py

@@ -1,17 +1,18 @@
-# Copyright (c) 2020-2024, Manfred Moitzi
+# Copyright (c) 2020-2025, Manfred Moitzi
 # License: MIT License
 from __future__ import annotations
 from typing import Sequence, Iterable, Optional, Iterator
 from enum import IntEnum
 import math
+import numpy as np
+
 from ezdxf.math import (
     Vec3,
     Vec2,
     Matrix44,
     X_AXIS,
     Y_AXIS,
     Z_AXIS,
-    AnyVec,
     UVec,
 )
 
@@ -38,6 +39,7 @@
     "bending_angle",
     "split_polygon_by_plane",
     "is_face_normal_pointing_outwards",
+    "is_vertex_order_ccw_3d",
 ]
 PI2 = math.pi / 2.0
 
@@ -82,7 +84,7 @@ def subdivide_face(
     """
     if len(face) < 3:
         raise ValueError("3 or more vertices required.")
-    
+
     len_face: int = len(face)
     mid_pos = Vec3.sum(face) / len_face
     subdiv_location: list[Vec3] = [
@@ -91,9 +93,7 @@ def subdivide_face(
 
     for index, vertex in enumerate(face):
         if quads:
-            yield vertex, subdiv_location[index], mid_pos, subdiv_location[
-                index - 1
-            ]
+            yield vertex, subdiv_location[index], mid_pos, subdiv_location[index - 1]
         else:
             yield subdiv_location[index - 1], vertex, mid_pos
             yield vertex, subdiv_location[index], mid_pos
@@ -378,9 +378,7 @@ def intersect_ray(self, origin: Vec3, direction: Vec3) -> Optional[Vec3]:
         """
         n = self.normal
         try:
-            weight = (self.distance_from_origin - n.dot(origin)) / n.dot(
-                direction
-            )
+            weight = (self.distance_from_origin - n.dot(origin)) / n.dot(direction)
         except ZeroDivisionError:
             return None
         return origin + (direction * weight)
@@ -646,9 +644,9 @@ def has_matrix_3d_stretching(m: Matrix44) -> bool:
     ux_mag_sqr = m.transform_direction(X_AXIS).magnitude_square
     uy = m.transform_direction(Y_AXIS)
     uz = m.transform_direction(Z_AXIS)
-    return not math.isclose(
-        ux_mag_sqr, uy.magnitude_square
-    ) or not math.isclose(ux_mag_sqr, uz.magnitude_square)
+    return not math.isclose(ux_mag_sqr, uy.magnitude_square) or not math.isclose(
+        ux_mag_sqr, uz.magnitude_square
+    )
 
 
 def spherical_envelope(points: Sequence[UVec]) -> tuple[Vec3, float]:
@@ -665,9 +663,7 @@ def spherical_envelope(points: Sequence[UVec]) -> tuple[Vec3, float]:
     return centroid, radius
 
 
-def inscribe_circle_tangent_length(
-    dir1: Vec3, dir2: Vec3, radius: float
-) -> float:
+def inscribe_circle_tangent_length(dir1: Vec3, dir2: Vec3, radius: float) -> float:
     """Returns the tangent length of an inscribe-circle of the given `radius`.
     The direction `dir1` and `dir2` define two intersection tangents,
     The tangent length is the distance from the intersection point of the
@@ -698,10 +694,10 @@ def bending_angle(dir1: Vec3, dir2: Vec3, normal=Z_AXIS) -> float:
 
 
 def any_vertex_inside_face(vertices: Sequence[Vec3]) -> Vec3:
-    """Returns a vertex from the "inside" of  the given face.
-    """
+    """Returns a vertex from the "inside" of  the given face."""
     # Triangulation is for concave shapes important!
     from ezdxf.math.triangulation import mapbox_earcut_3d
+
     it = mapbox_earcut_3d(vertices)
     return Vec3.sum(next(it)) / 3.0
 
@@ -718,12 +714,11 @@ def front_faces_intersect_face_normal(
     A counter-clockwise vertex order is assumed!
 
     """
+
     def is_face_in_front_of_detector(vertices: Sequence[Vec3]) -> bool:
         if len(vertices) < 3:
             return False
-        return any(
-            detector_plane.signed_distance_to(v) > abs_tol for v in vertices
-        )
+        return any(detector_plane.signed_distance_to(v) > abs_tol for v in vertices)
 
     # face-normal for counter-clockwise vertex order
     face_normal = safe_normal_vector(face)
@@ -764,6 +759,46 @@ def is_face_normal_pointing_outwards(
     This function does not check if the `faces` are a closed surface.
 
     """
-    return (
-        front_faces_intersect_face_normal(faces, face, abs_tol=abs_tol) % 2 == 0
-    )
+    return front_faces_intersect_face_normal(faces, face, abs_tol=abs_tol) % 2 == 0
+
+
+def is_vertex_order_ccw_3d(vertices: list[Vec3], normal: Vec3) -> bool:
+    """Returns ``True`` when the given 3D vertices have a counter-clockwise order around
+    the given normal vector.
+
+    Works for convex and concave shapes. Does not check or care if all vertices are
+    located in a flat plane or if the normal vector is really perpendicular to the
+    shape, but the result may be incorrect in that cases.
+
+    Args:
+        vertices (list): corner vertices of a flat shape (polygon)
+        normal (Vec3): normal vector of the shape
+
+    Raises:
+        ValueError: input has less than 3 vertices
+
+    """
+    if len(vertices) < 3:
+        raise ValueError("3 or more vertices required")
+
+    def signed_area() -> float:
+        # using the shoelace formula.
+        polygon = np.array(vertices)
+        if dom == 0:  # dominant X, use YZ plane
+            x, y = polygon[:, 1], polygon[:, 2]
+        elif dom == 1:  # dominant Y, use XZ plane
+            x, y = polygon[:, 0], polygon[:, 2]
+        else:  # dominant Z, use XY plane
+            x, y = polygon[:, 0], polygon[:, 1]
+        # returns twice the area, but only the sign is relevant for this use case
+        return np.dot(x, np.roll(y, -1)) - np.dot(y, np.roll(x, -1))
+
+    # The polygon is maybe concave, direct cross-product checks between
+    # adjacent edges are unreliable. Instead, use a projected signed area method.
+    # Find dominant axis:
+    abs_axis = abs(normal.x), abs(normal.y), abs(normal.z)
+    dom = abs_axis.index(max(abs_axis))
+
+    ccw = signed_area() > 0
+    dom_axis = normal[dom]
+    return dom_axis > 0 if ccw else dom_axis < 0

tests/test_06_math/test_614_construct_3d.py

@@ -1,7 +1,7 @@
 # Copyright (c) 2020, Manfred Moitzi
 # License: MIT License
 import pytest
-
+import math
 from ezdxf.math import (
     is_planar_face,
     Vec3,
@@ -21,15 +21,14 @@
     Matrix44,
     BarycentricCoordinates,
     linear_vertex_spacing,
+    is_vertex_order_ccw_3d,
 )
 
 from ezdxf.render.forms import square, circle
 
 REGULAR_FACE = Vec3.list([(0, 0, 0), (1, 0, 1), (1, 1, 1), (0, 1, 0)])
 IRREGULAR_FACE = Vec3.list([(0, 0, 0), (1, 0, 1), (1, 1, 0), (0, 1, 0)])
-REGULAR_FACE_WRONG_ORDER = Vec3.list(
-    [(0, 0, 0), (1, 1, 1), (1, 0, 1), (0, 1, 0)]
-)
+REGULAR_FACE_WRONG_ORDER = Vec3.list([(0, 0, 0), (1, 1, 1), (1, 0, 1), (0, 1, 0)])
 ONLY_COLINEAR_EDGES = Vec3.list([(0, 0, 0), (1, 0, 0), (2, 0, 0), (3, 0, 0)])
 REGULAR_FACE_WITH_COLINEAR_EDGE = Vec3.list(
     [(0, 0, 0), (1, 0, 0), (2, 0, 0), (3, 0, 0), (1.5, 2.0, 0)]
@@ -115,9 +114,7 @@ def test_intersecting_rays_return_one_tuple(self, ray1, ray2):
         assert bool(result) is True
         assert result == (Vec3(0, 0, 0),)
 
-    def test_not_intersecting_and_not_parallel_rays_return_two_tuple(
-        self, ray1, ray2
-    ):
+    def test_not_intersecting_and_not_parallel_rays_return_two_tuple(self, ray1, ray2):
         line3 = (Vec3(0, 0, 1), Vec3(0, 1, 1))
         result = intersection_ray_ray_3d(ray1, line3)
         assert len(result) == 2
@@ -161,22 +158,16 @@ def line4(self):
         return Vec3(2, -1, 0), Vec3(2, 1, 0)
 
     def test_real_intersecting_lines(self, line1, line2):
-        assert intersection_line_line_3d(line1, line2, virtual=False).isclose(
-            (1, 0, 0)
-        )
+        assert intersection_line_line_3d(line1, line2, virtual=False).isclose((1, 0, 0))
 
     def test_virtual_intersecting_lines(self, line1, line3):
-        assert intersection_line_line_3d(line1, line3, virtual=True).isclose(
-            (3, 0, 0)
-        )
+        assert intersection_line_line_3d(line1, line3, virtual=True).isclose((3, 0, 0))
 
     def test_not_intersecting_lines(self, line1, line3):
         assert intersection_line_line_3d(line1, line3, virtual=False) is None
 
     def test_touching_lines_do_intersect(self, line1, line4):
-        assert intersection_line_line_3d(line1, line4, virtual=False).isclose(
-            (2, 0, 0)
-        )
+        assert intersection_line_line_3d(line1, line4, virtual=False).isclose((2, 0, 0))
 
     @pytest.mark.parametrize(
         "p2", [(4, 0), (0, 4), (4, 4)], ids=["horiz", "vert", "diag"]
@@ -336,5 +327,56 @@ def test_correct_spacing_in_Q3(self, count):
             assert vertices[x].isclose((-x, -x, -x))
 
 
+I_BEAM = Vec3.list(
+    [
+        (0, 0),
+        (3, 0),
+        (3, 1),
+        (2, 1),
+        (2, 2),
+        (3, 2),
+        (3, 3),
+        (0, 3),
+        (0, 2),
+        (1, 2),
+        (1, 1),
+        (0, 1),
+    ]
+)
+
+
+class TestIsVertexOrderCCW:
+    def test_xy_plane(self):
+        assert is_vertex_order_ccw_3d(I_BEAM, Vec3(0, 0, 1)) is True
+
+    def test_xy_plane_inv(self):
+        assert is_vertex_order_ccw_3d(I_BEAM, Vec3(0, 0, -1)) is False
+
+    def test_yz_plane_up(self):
+        m = Matrix44.x_rotate(math.pi / 2)
+        vertices = list(m.transform_vertices(I_BEAM))
+        assert is_vertex_order_ccw_3d(vertices, Vec3(0, 1, 0)) is True
+
+    def test_yz_plane_inv(self):
+        m = Matrix44.x_rotate(math.pi / 2)
+        vertices = list(m.transform_vertices(I_BEAM))
+        assert is_vertex_order_ccw_3d(vertices, Vec3(0, -1, 0)) is False
+
+    def test_xz_plane_up(self):
+        m = Matrix44.y_rotate(math.pi / 2)
+        vertices = list(m.transform_vertices(I_BEAM))
+        assert is_vertex_order_ccw_3d(vertices, Vec3(1, 0, 0)) is True
+
+    def test_xz_plane_inv(self):
+        m = Matrix44.y_rotate(math.pi / 2)
+        vertices = list(m.transform_vertices(I_BEAM))
+        assert is_vertex_order_ccw_3d(vertices, Vec3(-1, 0, 0)) is False
+
+    def test_square_xy_plane(self):
+        square = Vec3.list([(0, 0), (1, 0), (1, 1), (0, 1)])
+        assert is_vertex_order_ccw_3d(square, Vec3(0, 0, 1)) is True
+        assert is_vertex_order_ccw_3d(square, Vec3(0, 0, -1)) is False
+
+
 if __name__ == "__main__":
     pytest.main([__file__])