GouMinghao / Geometry3D

Geometry3D: 3D Computational Geometrics Library
GNU General Public License v3.0
88 stars 23 forks source link

Normals for polys where first 3 points are in a line #17

Open FriendOfThePigeon opened 1 year ago

FriendOfThePigeon commented 1 year ago

When you create an intersection, it seems that sometimes the ConvexPolygons that are created have 3 points along the first side. When that happens, the creation of the normal fails with "float division by zero". This is maybe not the most elegant solution but, since the polygon is convex, there must be another set of 3 points that are not in a line.

I can provide example code that shows the problem, but it doesn't occur every time, so creating a test is difficult.

FriendOfThePigeon commented 1 year ago

Code to replicate:

#!/usr/bin/env python

import sys
from itertools import chain

from Geometry3D import Vector, ConvexPolygon, ConvexPolyhedron, Point, intersection

def tri(p1, p2, p3):
    return ConvexPolygon((p1, p2, p3))

def quad(p1, p2, p3, p4):
    yield tri(p1, p2, p3)
    yield tri(p3, p4, p1)

def cph_opposite_faces(face1, face2):
    (p1, p2, p3, p4, p5, p6, p7, p8) = [face[i] for face in [face1, face2] for i in range(4)]
    return ConvexPolyhedron(tuple(chain(
        quad(p4, p3, p2, p1),
        quad(p1, p5, p8, p4),
        quad(p5, p1, p2, p6),
        quad(p5, p6, p7, p8),
        quad(p4, p8, p7, p3),
        quad(p2, p3, p7, p6) )))

def opposite_faces(cls, name, face1, face2, material=None):
    return Form(
        name, 
        _cph_opposite_faces(face1, face2),
        material=material)

def axis_aligned(pos, size):
    # An axis-aligned cuboid, defined by position and size
    assert all(size[i] > 0 for i in range(3)), 'Size must be > 0 in all dims'
    (x1, x2, y1, y2, z1, z2) = list(chain.from_iterable([pos[i], pos[i] + size[i]] for i in range(3)))
    result = cph_opposite_faces(
            (Point(x1, y1, z1), Point(x1, y1, z2), Point(x2, y1, z2), Point(x2, y1, z1)),
            (Point(x1, y2, z1), Point(x1, y2, z2), Point(x2, y2, z2), Point(x2, y2, z1))
            )
    return result

one = axis_aligned((50, 50, -5000), (10000, 10000, 10000))
two = axis_aligned((0, 0, 0), (100, 100, 100))

result = intersection(one, two)  # Sometimes this fails, with "float division by zero"
print(str(result))