Skip to content

Commit

Permalink
Fix flow_to_exit() in more cases
Browse files Browse the repository at this point in the history 
fixes flatsurf#259
  • Loading branch information
@saraedum
saraedum committed Oct 10, 2023
1 parent 96376a1 commit 7547606
Show file tree
Hide file tree
Showing 2 changed files with 79 additions and 68 deletions.
109 changes: 41 additions & 68 deletions flatsurf/geometry/categories/euclidean_polygons.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1564,7 +1564,7 @@ def contains_point(self, point, translation=None):
point, translation=translation
).is_inside()

def flow_to_exit(self, point, direction):
def flow_to_exit(self, point, direction, stop_at_vertex=False):
r"""
Flow a point in the direction of holonomy until the point leaves the
polygon. Note that ValueErrors may be thrown if the point is not in the
Expand All @@ -1590,77 +1590,50 @@ def flow_to_exit(self, point, direction):
sage: P.flow_to_exit(vector((2, 1)), vector((0, 1)))
((2, 3), point positioned on vertex 6 of polygon)
"""
# TODO: Maybe this should be rewritten to be less cryptic
# and work more robustly for non-strictly convex polygons.
from flatsurf.geometry.polygon import PolygonPosition
if not direction:
raise ValueError("direction must be non-zero")

V = self.base_ring().fraction_field() ** 2
if direction == V.zero():
raise ValueError("Zero vector provided as direction.")
v0 = self.vertex(0)
for i in range(len(self.vertices())):
e = self.edge(i)
from sage.all import matrix
vertices = self.vertices()

m = matrix([[e[0], -direction[0]], [e[1], -direction[1]]])
try:
ret = m.inverse() * (point - v0)
s = ret[0]
t = ret[1]
# What if the matrix is non-invertible?
first_intersection = None

# Answer: You'll get a ZeroDivisionError which means that the edge is parallel
# to the direction.
for v in range(len(vertices)):
segment = vertices[v], vertices[(v + 1) % len(vertices)]

# s is location it intersects on edge, t is the portion of the direction to reach this intersection
if t > 0 and 0 <= s and s <= 1:
# The ray passes through edge i.
if s == 1:
# exits through vertex i+1
v0 = v0 + e
return v0, PolygonPosition(
PolygonPosition.VERTEX,
vertex=(i + 1) % len(self.vertices()),
)
if s == 0:
# exits through vertex i
return v0, PolygonPosition(
PolygonPosition.VERTEX, vertex=i
)
# exits through vertex i
# exits through interior of edge i
prod = t * direction
return point + prod, PolygonPosition(
PolygonPosition.EDGE_INTERIOR, edge=i
)
except ZeroDivisionError:
# Here we know the edge and the direction are
# parallel or anti-parallel.
from flatsurf.geometry.euclidean import is_parallel, is_anti_parallel
if point == v0 or (is_parallel(e, point - v0) and is_anti_parallel(e, point - v0 + e)):
# In this case point lies on the edge.
# We need to work out which direction to move in.

if (point - v0).is_zero() or is_parallel(e, point - v0):
# exits through vertex i+1
return self.vertex(i + 1), PolygonPosition(
PolygonPosition.VERTEX,
vertex=(i + 1) % len(self.vertices()),
)
else:
# exits through vertex i
return v0, PolygonPosition(
PolygonPosition.VERTEX, vertex=i
)
pass
v0 = v0 + e
# Our loop has terminated. This can mean one of several errors...
pos = self.get_point_position(point)
if pos.is_outside():
raise ValueError("Started with point outside polygon")
raise ValueError(
"Point on boundary of polygon and direction not pointed into the polygon."
)
from flatsurf.geometry.euclidean import ray_segment_intersection
intersection = ray_segment_intersection(point, direction, segment)

if intersection is None:
continue

if isinstance(intersection, tuple):
if intersection[0] != point:
# The flow overlaps with this edge but it hits
# a vertex before it gets here.
continue

intersection = intersection[1]
assert intersection != point

if intersection == point:
continue

# TODO: Implement stop_at_vertex == False. (But what
# should be the default actually? I think the old
# version was just buggy so there is no legacy to
# follow.)
from flatsurf.geometry.euclidean import time_on_ray
if first_intersection is None or time_on_ray(point, direction, first_intersection) > time_on_ray(point, direction, intersection):
first_intersection = intersection

if first_intersection is not None:
# TODO: Do not call expensive get_point_position.
return first_intersection, self.get_point_position(first_intersection)

if self.get_point_position(point).is_outside():
raise ValueError("Cannot flow from point outside of polygon")

raise ValueError("Cannot flow from point on boundary if direction points out of the polygon")

def flow_map(self, direction):
r"""
Expand Down
38 changes: 38 additions & 0 deletions flatsurf/geometry/euclidean.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -383,6 +383,44 @@ def line_intersection(p1, p2, q1, q2):
return p1 + (m.inverse() * (q1 - p1))[0] * (p2 - p1)


def time_on_ray(p, direction, q):
if direction[0]:
return (q[0] - p[0]) / direction[0]
else:
return (q[1] - p[1]) / direction[1]


def ray_segment_intersection(p, direction, segment):
intersection = line_intersection(p, p + direction, *segment)

if intersection is None:
# ray and segment are parallel.
t0 = time_on_ray(p, direction, segment[0])
t1 = time_on_ray(p, direction, segment[1])

if t1 < t0:
t0, t1 = t1, t0

if t1 < 0:
return None

if t1 == 0:
return p

if t0 < 0:
return (p, p + t1 * direction)

return (p + t0 * direction, p + t1 * direction)

if time_on_ray(p, direction, intersection) < 0:
return None

if ccw(segment[0] - p, direction) * ccw(segment[1] - p, direction) == 1:
return None

return intersection


def is_segment_intersecting(e1, e2):
r"""
Return whether the segments ``e1`` and ``e2`` intersect.
Expand Down

0 comments on commit 7547606

Please sign in to comment.