Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix some issues with Sketch.hull() #1344

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Fix some issues with Sketch.hull() #1344

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
c84b682
Fix hull
TerranceN Jun 10, 2023
c519f99
Get tests to pass by fixing another edge case
TerranceN Jun 10, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
71 changes: 49 additions & 22 deletions cadquery/hull.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,8 +70,8 @@ def __init__(self, c: Point, r: float, a1: float, a2: float):
self.a1 = a1
self.a2 = a2

self.s = Point(r * cos(a1), r * sin(a1))
self.e = Point(r * cos(a2), r * sin(a2))
self.s = Point(c.x + r * cos(a1), c.y + r * sin(a1))
self.e = Point(c.x + r * cos(a2), c.y + r * sin(a2))
self.ac = 2 * pi - (a1 - a2)


Expand Down Expand Up @@ -104,10 +104,12 @@ def convert_and_validate(edges: Iterable[Edge]) -> Tuple[List[Arc], List[Point]]
r = e.radius()
a1, a2 = e._bounds()

arcs.add(Arc(Point(c.x, c.y), r, a1, a2))
arc = Arc(Point(c.x, c.y), r, a1, a2)
arcs.add(arc)
points.update((arc.s, arc.e))

else:
raise ValueError("Unsupported geometry {gt}")
raise ValueError(f"Unsupported geometry {gt}")

return list(arcs), list(points)

Expand Down Expand Up @@ -187,6 +189,9 @@ def _pt_arc(p: Point, a: Arc) -> Tuple[float, float, float, float]:
xc, yc = a.c.x, a.c.y
dx, dy = x - xc, y - yc
l = sqrt(dx ** 2 + dy ** 2)
ld = l ** 2 - r ** 2
if (ld < 0.001):
return x, y, x, y
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is when the point is on the arc's circle, so I return a no-length segment (which is later ignored). The point might not necessarily be on the arc though, so this should return a segment to the start or end of the arc.

This is taken care of elsewhere because I add the start and end points of the arc to the list of points in convert_and_validate, so if the next segment is between the point and a start or end point of an arc, it's handled elsewhere as a pt_pt comparison.


x1 = r ** 2 / l ** 2 * dx - r / l ** 2 * sqrt(l ** 2 - r ** 2) * dy + xc
y1 = r ** 2 / l ** 2 * dy + r / l ** 2 * sqrt(l ** 2 - r ** 2) * dx + yc
Expand All @@ -195,33 +200,38 @@ def _pt_arc(p: Point, a: Arc) -> Tuple[float, float, float, float]:

return x1, y1, x2, y2


def pt_arc(p: Point, a: Arc) -> Tuple[float, Segment]:
def pt_arc(angle: float, p: Point, a: Arc) -> Tuple[float, Segment]:

x, y = p.x, p.y
x1, y1, x2, y2 = _pt_arc(p, a)

if ((x1, y1) == (x2, y2)):
return inf, Segment(Point(inf, inf), Point(inf, inf))

angles = atan2p(x1 - x, y1 - y), atan2p(x2 - x, y2 - y)
points = Point(x1, y1), Point(x2, y2)
ix = int(argmin(angles))
ix = int(argmin([a + 2*pi if a < angle else a for a in angles]))

return angles[ix], Segment(p, points[ix])


def arc_pt(a: Arc, p: Point) -> Tuple[float, Segment]:
def arc_pt(angle: float, a: Arc, p: Point) -> Tuple[float, Segment]:

x, y = p.x, p.y
x1, y1, x2, y2 = _pt_arc(p, a)

if ((x1, y1) == (x2, y2)):
return inf, Segment(Point(inf, inf), Point(inf, inf))

angles = atan2p(x - x1, y - y1), atan2p(x - x2, y - y2)
points = Point(x1, y1), Point(x2, y2)

ix = int(argmax(angles))
ix = int(argmax([a + 2*pi if a < angle else a for a in angles]))
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When the angle loops around back past 0 radians, this caused the argmin/argmax of pt_arc and arc_pt to choose the wrong side of the arc to continue from.


return angles[ix], Segment(points[ix], p)


def arc_arc(a1: Arc, a2: Arc) -> Tuple[float, Segment]:
def arc_arc(angle: float, a1: Arc, a2: Arc) -> Tuple[float, Segment]:

r1 = a1.r
xc1, yc1 = a1.c.x, a1.c.y
Expand Down Expand Up @@ -292,24 +302,36 @@ def arc_arc(a1: Arc, a2: Arc) -> Tuple[float, Segment]:
Segment(Point(x12, y12), Point(x22, y22)),
)

return angles[ix], segments[ix]
segment_angle = atan2p(segments[ix].b.x - xc2, segments[ix].b.y - yc2)
if (segment_angle < a2.a1 or segment_angle > a2.a2):
return inf, Segment(Point(inf, inf), Point(inf, inf))
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand this function that well, so this might have been from something else wrong and is now redundant, but somehow it was selecting a point not actually on the arc, so this if statement throws that out if that happens.

I'd bet this isn't the best way to handle it and there's more edge cases that would fail here.


return angles[ix], segments[ix]

def get_angle(current: Entity, e: Entity) -> Tuple[float, Segment]:

if current is e:
def get_angle(angle: float, current: Entity, e: Entity) -> Tuple[float, Segment]:
def center(e: Entity):
if isinstance(e, Point):
return e
return e.c
if center(current) == center(e):
if (isinstance(current, Arc) and isinstance(e, Arc)):
if (current.s == e.e and not e.s == current.e):
return pt_pt(current.e, e.s)
if (e.s == current.e and not current.s == e.e):
return (current.a2 + pi) % (2 * pi), Segment(Point(inf, inf), Point(inf, inf))
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These changes allow a circle that's broken into multiple arcs continue around the circle

return inf, Segment(Point(inf, inf), Point(inf, inf))

if isinstance(current, Point):
if isinstance(e, Point):
return pt_pt(current, e)
else:
return pt_arc(current, e)
return pt_arc(angle, current, e)
else:
if isinstance(e, Point):
return arc_pt(current, e)
return arc_pt(angle, current, e)
else:
return arc_arc(current, e)
return arc_arc(angle, current, e)


def update_hull(
Expand Down Expand Up @@ -339,7 +361,8 @@ def finalize_hull(hull: Hull) -> Wire:
for el_p, el, el_n in zip(hull, hull[1:], hull[2:]):

if isinstance(el, Segment):
rv.append(Edge.makeLine(Vector(el.a.x, el.a.y), Vector(el.b.x, el.b.y)))
if (not inf in [el.a.x, el.a.y, el.b.x, el.b.y]):
rv.append(Edge.makeLine(Vector(el.a.x, el.a.y), Vector(el.b.x, el.b.y)))
elif (
isinstance(el, Arc)
and isinstance(el_p, Segment)
Expand All @@ -348,9 +371,10 @@ def finalize_hull(hull: Hull) -> Wire:
a1 = degrees(atan2p(el_p.b.x - el.c.x, el_p.b.y - el.c.y))
a2 = degrees(atan2p(el_n.a.x - el.c.x, el_n.a.y - el.c.y))

rv.append(
Edge.makeCircle(el.r, Vector(el.c.x, el.c.y), angle1=a1, angle2=a2)
)
if (a1 != a2):
rv.append(
Edge.makeCircle(el.r, Vector(el.c.x, el.c.y), angle1=a1, angle2=a2)
)

el1 = hull[1]
if isinstance(el, Segment) and isinstance(el_n, Arc) and isinstance(el1, Segment):
Expand Down Expand Up @@ -385,14 +409,17 @@ def find_hull(edges: Iterable[Edge]) -> Wire:
current_angle = 0.0
finished = False

iterations = 0

# march around
while not finished:
while not finished and iterations < len(entities):
iterations += 1
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

While working on this it was very easy for this to loop infinitely. Setting a sensible maximum number of iterations (if we've iterated more than the number of entities, something has gone wrong) seems better than having the system hang in an infinite loop.


angles = []
segments = []

for e in entities:
angle, segment = get_angle(current_e, e)
angle, segment = get_angle(current_angle, current_e, e)
angles.append(angle if angle >= current_angle else inf)
segments.append(segment)

Expand Down
2 changes: 1 addition & 1 deletion tests/test_hull.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ def test_hull():

h = hull.find_hull(edges)

assert len(h.Vertices()) == 11
assert len(h.Vertices()) == 10
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Somehow my changes made this hull have one fewer vertices. Based on how the points and arcs are laid out, I'm guessing there were multiple possible next entities with the same angle and my hull changes made it skip one.

There should probably be a 'tie-breaker' for choosing the smallest angle so that the system is more deterministic. (either always choosing the closer to ensure the loop always ends up closed, or choosing the further, but we'd have to be more careful about choosing the starting entity, lest we start on an entity we later skip and fail to close the loop)

assert h.IsClosed()
assert h.isValid()

Expand Down