Clipping!

src/geometry_utils/intersections.jl

@@ -220,4 +220,68 @@ function intersection_of_ray_with_edge(p, q, a, b)
     else
         return (p1 + t * (q1 - p1), p2 + t * (q2 - p2))
     end
-end
+end
+
+"""
+    liang_barsky(a, b, c, d, p, q)
+
+Applies the Liang-Barsky algorithm to find the intersection of the line segment from `p` to `q`
+with the rectangle from `[a, b] × [c, d]`. Returns `(u, v)`, where:
+
+- If there is an intersection, then `u` and `v` are the coordinates of the two intersections. 
+- If there is no intersection, then `u = v = (NaN, NaN)`.
+"""
+function liang_barsky(a, b, c, d, p, q)
+    t1 = 0.0
+    t2 = 1.0
+    px, py = _getxy(p)
+    qx, qy = _getxy(q)
+    Δx = qx - px
+    t1, t2, inside = liang_barsky_clipper(-Δx, px - a, t1, t2)
+    if inside
+        t1, t2, inside = liang_barsky_clipper(Δx, b - px, t1, t2)
+        if inside
+            Δy = qy - py
+            t1, t2, inside = liang_barsky_clipper(-Δy, py - c, t1, t2)
+            if inside
+                t1, t2, inside = liang_barsky_clipper(Δy, d - py, t1, t2)
+                if inside
+                    if t2 < 1
+                        qx = px + t2 * Δx
+                        qy = py + t2 * Δy
+                    end
+                    if t1 > 0
+                        px = px + t1 * Δx
+                        py = py + t1 * Δy
+                    end
+                end
+            end
+        end
+    end
+    if inside
+        return (px, py), (qx, qy)
+    else
+        return (NaN, NaN), (NaN, NaN)
+    end
+end
+function liang_barsky_clipper(p, q, t1, t2)
+    inside = true
+    if p < 0
+        r = q / p
+        if r > t2
+            inside = false
+        elseif r > t1
+            t1 = r
+        end
+    elseif p > 0
+        r = q / p
+        if r < t1
+            inside = false
+        elseif r < t2
+            t2 = r
+        end
+    elseif q < 0
+        inside = false
+    end
+    return t1, t2, inside
+end

src/geometry_utils/sutherland_hodgman.jl

@@ -1,7 +1,7 @@
 struct Polygon{T,V,P} <: AbstractVector{T}
     vertices::V
     points::P
-    function Polygon(vertices::V, points::P) where {V,P}
+    @inline function Polygon(vertices::V, points::P) where {V,P}
         p = get_point(points, vertices[begin])
         T = typeof(p)
         if vertices[begin] ≠ vertices[end]

src/voronoi/coordinates.jl

@@ -70,25 +70,26 @@ end
 """
     _get_ray(vorn, i, boundary_index)
 
-Extracts the ray from the `i`th polygon of `vorn` corresponding to the `boundary_index`th boundary index. 
+Extracts the ray from the `i`th polygon of `vorn` corresponding to the `boundary_index`, where `boundary_index` 
+here means that `get_polygon(vorn, i)[boundary_index]` is a boundary index.
 The returned points are given in the form `(p, q)`, defining the oriented line `pq` such that the line 
 is in the direction of infinity.
 """
 function _get_ray(vorn, i, boundary_index)
-    ghost_tri = get_circumcenter_to_triangle(vorn, boundary_index)
+    C = get_polygon(vorn, i)
+    ghost_tri = get_circumcenter_to_triangle(vorn, C[boundary_index])
     u, v, _ = indices(ghost_tri) # w is the ghost vertex
     p, q = get_generator(vorn, u, v)
     px, py = _getxy(p)
     qx, qy = _getxy(q)
     mx, my = (px + qx) / 2, (py + qy) / 2
     m = (mx, my)
-    C = get_polygon(vorn, i)
-    is_first = is_first_boundary_index(C, i)
+    is_first = is_first_boundary_index(C, boundary_index)
     if is_first
-        prev_index = previndex_circular(C, i)
+        prev_index = previndex_circular(C, boundary_index)
         r = get_polygon_point(vorn, C[prev_index])
     else
-        next_index = nextindex_circular(C, i)
+        next_index = nextindex_circular(C, boundary_index)
         r = get_polygon_point(vorn, C[next_index])
     end
     if r == m # It's possible for the circumcenter to lie on the edge and exactly at the midpoint (e.g. [(0.0,1.0),(-1.0,2.0),(-2.0,-1.0)]). In this case, just rotate 
@@ -107,6 +108,77 @@ function _get_ray(vorn, i, boundary_index)
     return m, r
 end
 
+"""
+    grow_polygon_outside_of_box(vorn::VoronoiTessellation, i, bounding_box)
+
+Truncates unbounded edges of the `i`th polygon of `vorn`, assumed to be unbounded,
+so that the line connecting the truncated unbounded edges is entirely outside 
+of the polygon. The method of growth is iterative, utilising the Liang-Barsky algorithm 
+at each stage while we translate the line. The returned polygon does not satisfy 
+`P[begin] == P[end]`.
+"""
+function grow_polygon_outside_of_box(vorn::VoronoiTessellation, i, bounding_box)
+    a, b, c, d = bounding_box
+    vertices = get_polygon(vorn, i)
+    new_vertices, new_points, boundary_indices = get_new_polygon_indices(vorn, vertices)
+    inside = true
+    t = 1.0 # don't do 0.5 so we get t = 1 later, else we get duplicated vertices for polygons completely outside of the box
+    u, v = boundary_indices
+    u_m, u_r = _get_ray(vorn, i, u)
+    v_m, v_r = _get_ray(vorn, i, v)
+    u_mx, u_my = _getxy(u_m)
+    u_rx, u_ry = _getxy(u_r)
+    v_mx, v_my = _getxy(v_m)
+    v_rx, v_ry = _getxy(v_r)
+    p = (0.0, 0.0)
+    q = (0.0, 0.0)
+    while inside
+        t *= 2.0
+        p = (u_mx + t * (u_rx - u_mx), u_my + t * (u_ry - u_my))
+        q = (v_mx + t * (v_rx - v_mx), v_my + t * (v_ry - v_my))
+        int1, int2 = liang_barsky(a, b, c, d, p, q)
+        outside = all(isnan, int1) && all(isnan, int2)
+        inside = !outside
+    end
+    new_points[u] = p
+    new_points[v] = q
+    new_vertices[u] = u
+    new_vertices[v] = v
+    return new_vertices, new_points
+end
+
+"""
+    get_new_polygon_indices(vorn, vertices)
+
+Given an unbounded Voronoi polygon from `vorn` with `vertices`, returns 
+`(new_vertices, new_points, boundary_indices)`, where `new_vertices` are vertices 
+mapping to the points in `new_points`, which is just a vector of all the polygon points, 
+and `boundary_indices` is a `Tuple` of the indices of the points in `new_points` that correspond 
+to points out at infinity.
+"""
+function get_new_polygon_indices(vorn, vertices)
+    new_points = NTuple{2,Float64}[]
+    sizehint!(new_points, length(vertices))
+    new_vertices = similar(vertices, length(vertices) - 1)
+    boundary_indices = (0, 0)
+    for i in firstindex(vertices):(lastindex(vertices)-1)
+        v = vertices[i]
+        if is_boundary_index(v)
+            is_first = is_first_boundary_index(vertices, i)
+            if is_first
+                boundary_indices = (i, boundary_indices[2])
+            else
+                boundary_indices = (boundary_indices[1], i)
+            end
+            push!(new_points, (NaN, NaN))
+            new_vertices[i] = v
+        else
+            push!(new_points, _getxy(get_polygon_point(vorn, v)))
+            new_vertices[i] = length(new_points)
+        end
+    end
+    return new_vertices, new_points, boundary_indices
+end
 
 """
     get_bounded_polygon_coordinates(vorn::VoronoiTessellation, i, bounding_box)
@@ -144,3 +216,19 @@ function get_unbounded_polygon_coordinates(vorn::VoronoiTessellation, i, boundin
         return clip_unbounded_polygon_to_bounding_box(vorn, i, bounding_box)
     end
 end
+
+"""
+    clip_unbounded_polygon_to_bounding_box(vorn::VoronoiTessellation, i, bounding_box)
+
+Clips the `i`th polygon of `vorn` to the bounding box, assuming it is unbounded.
+The unbounded polygon is truncated so that the line connecting the unbounded edges 
+is outside of the bounded box, and then the Sutherland-Hodgman algorithm 
+is used to clip the resulting polygon to the bounding box.
+"""
+function clip_unbounded_polygon_to_bounding_box(vorn::VoronoiTessellation, i, bounding_box)
+    new_vertices, new_points = grow_polygon_outside_of_box(vorn, i, bounding_box)
+    clip_vertices = (1, 2, 3, 4)
+    a, b, c, d = bounding_box
+    clip_points = ((a, c), (b, c), (b, d), (a, d))
+    return clip_polygon(new_vertices, new_points, clip_vertices, clip_points)
+end

test/geo_utils.jl

@@ -1094,4 +1094,42 @@ end
       b = (-2.0, 0.0)
       r = DT.intersection_of_ray_with_edge(p, q, a, b)
       @test collect(r) ≈ [0.0, 0.0]
-end
+end
+
+@testset "Liang-Barsky algorithm" begin
+      p, q = (-7.0, 5.0), (1.0, 8.0)
+      a, b, c, d = 0.0, 8.0, 0.0, 4.0
+      u, v = DT.liang_barsky(a, b, c, d, p, q)
+      @test all(isnan, u) && all(isnan, v)
+      p, q = (-3.0, 1.0), (10.0, 5.0)
+      u, v = DT.liang_barsky(a, b, c, d, p, q)
+      @test collect(u) ≈ [0, 1.9230769230769]
+      @test collect(v) ≈ [6.75, 4.0]
+      p, q = (0.0, -2.0), (0.0, 6.0)
+      u, v = DT.liang_barsky(a, b, c, d, p, q)
+      @test collect(u) ≈ [0, 0]
+      @test collect(v) ≈ [0, 4]
+      p, q = (2.0, 6.0), (-2.0, 2.0)
+      u, v = DT.liang_barsky(a, b, c, d, p, q)
+      @test u == v && (collect(u) ≈ [0, 4])
+      p, q = (10.0, 6.0), (-2.0, 6.0)
+      u, v = DT.liang_barsky(a, b, c, d, p, q)
+      @test all(isnan, u) && all(isnan, v)
+      p, q = (4.0, 6.0), (4.0, -2.0)
+      u, v = DT.liang_barsky(a, b, c, d, p, q)
+      @test collect(u) ≈ [4.0, 4.0]
+      @test collect(v) ≈ [4.0, 0.0]
+      p, q = (2.0, 6.0), (10.0, -2.0)
+      u, v = DT.liang_barsky(a, b, c, d, p, q)
+      @test collect(u) ≈ [4.0, 4.0]
+      @test collect(v) ≈ [8.0, 0.0]
+      a, b, c, d = 4, 10, 2, 8
+      u, v = DT.liang_barsky(a, b, c, d, p, q)
+      @test collect(u) ≈ [4.0, 4.0]
+      @test collect(v) ≈ [6.0, 2.0]
+      p, q = (2.0, 6.0), (14.0, 8.0)
+      u, v = DT.liang_barsky(a, b, c, d, p, q)
+      @test collect(u) ≈ [4, 6 + 1 / 3]
+      @test collect(v) ≈ [10, 7 + 1 / 3]
+end
+