Skip to content

Commit

Permalink
Clipping Voronoi tessellations to arbitrary bounding boxes (#76)
Browse files Browse the repository at this point in the history
  • Loading branch information
@DanielVandH
DanielVandH authored Aug 10, 2023
1 parent dc18f8b commit 7ddb44a
Show file tree
Hide file tree
Showing 15 changed files with 920 additions and 204 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "DelaunayTriangulation"
uuid = "927a84f5-c5f4-47a5-9785-b46e178433df"
authors = ["Daniel VandenHeuvel <[email protected]>"]
version = "0.8.3"
version = "0.8.4"

[deps]
DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
Expand Down
Binary file modified docs/src/tessellations/figs/unbounded.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
3 changes: 3 additions & 0 deletions src/DelaunayTriangulation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -191,6 +191,7 @@ export convert_boundary_points_to_indices
include("geometry_utils/polygons.jl")
include("geometry_utils/polylabel.jl")
include("geometry_utils/intersections.jl")
include("geometry_utils/sutherland_hodgman.jl")

const polylabel = pole_of_inaccessibility

Expand All @@ -205,8 +206,10 @@ include("voronoi/main.jl")
include("voronoi/unbounded_construction.jl")
include("voronoi/clipped_construction.jl")
include("voronoi/lloyd.jl")
include("voronoi/coordinates.jl")

export voronoi
export centroidal_smooth
export get_polygon_coordinates

end
164 changes: 21 additions & 143 deletions src/data_structures/voronoi/voronoi.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,9 @@ while the values are the coordinates; we need to use a `Dict` in case the triang
The points defining the vertices of the polygons. The points are not guaranteed to be unique if a circumcenter
appears on the boundary and you are considering a clipped tessellation.
See also [`get_polygon_coordinates`](@ref).
- `polygons::Dict{I, Vector{I}}`
A `Dict` mapping a polygon index (same as a generator index) to the vertices of the polygon. The polygons are given in counter-clockwise order,
Expand Down Expand Up @@ -150,136 +153,6 @@ function get_triangle_to_circumcenter(vor::VoronoiTessellation, T)
end
end

"""
get_polygon_coordinates(vor::VoronoiTessellation, j, bbox=nothing)
Gets the coordinates of the `j`th polygon of `vor`. If `bbox` is given, then the coordinates are clipped to the bounding box `bbox`, with
`bbox` a `Tuple` of the form `(xmin, xmax, ymin, ymax)` so that the bounding box is `[xmin, xmax] × [ymin, ymax]`.
See also [`polygon_bounds`](@ref) for `bbox`.
"""
function get_polygon_coordinates(vorn::VoronoiTessellation, j, bbox=nothing)
C = get_polygon(vorn, j)
F = number_type(vorn)
coords = Vector{NTuple{2,F}}(undef, length(C) - 1)
unbounded_indices = (0, 0)
for i in firstindex(C):(lastindex(C)-1)
if !is_boundary_index(C[i])
coords[i] = get_polygon_point(vorn, C[i])
else
ghost_tri = get_circumcenter_to_triangle(vorn, C[i])
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)
@assert bbox[1] ≤ px ≤ bbox[2] && bbox[1] ≤ qx ≤ bbox[2] && bbox[3] ≤ py ≤ bbox[4] && bbox[3] ≤ qy ≤ bbox[4] "The bounding box is not large enough to contain the circumcenter."
m = (px + qx) / 2, (py + qy) / 2
is_first = is_first_boundary_index(C, i)
if is_first
prev_index = previndex_circular(C, i)
r = get_polygon_point(vorn, C[prev_index])
else
next_index = nextindex_circular(C, i)
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
mx, my = _getxy(m)
dx, dy = qx - mx, qy - my
rotated_dx, rotated_dy = -dy, dx
r = mx + rotated_dx, my + rotated_dy
if is_right(point_position_relative_to_line(p, q, r))
rotated_dx, rotated_dy = dy, -dx
r = mx + rotated_dx, my + rotated_dy
end
end
r = _getxy(r)
if is_left(point_position_relative_to_line(p, q, r))
intersection = intersection_of_ray_with_bounding_box(m, r, bbox[1], bbox[2], bbox[3], bbox[4])
else
intersection = intersection_of_ray_with_bounding_box(r, m, bbox[1], bbox[2], bbox[3], bbox[4])
end
coords[i] = intersection
if is_first
unbounded_indices = (i, unbounded_indices[2])
else
unbounded_indices = (unbounded_indices[1], i)
end
end
end
if all(≠(0), unbounded_indices) # Need to check if we need to insert the corner point of the polygon
xmin, xmax, ymin, ymax = bbox
c1 = coords[unbounded_indices[1]]
c2 = coords[unbounded_indices[2]]
side1 = identify_side(c1, xmin, xmax, ymin, ymax)
side2 = identify_side(c2, xmin, xmax, ymin, ymax)
if side1 ≠ side2 # There is a better way to do this if we treat sides as numbers 0, 1, 2, 3 rather than symbols, but it doesn't really matter.
if side1 == :bottom && side2 == :right
corner = xmax, ymin
insert!(coords, unbounded_indices[1] + 1, corner)
elseif side1 == :bottom && side2 == :top
corner1 = xmax, ymin
corner2 = xmax, ymax
insert!(coords, unbounded_indices[1] + 1, corner1)
insert!(coords, unbounded_indices[1] + 2, corner2)
elseif side1 == :bottom && side2 == :left
corner1 = xmax, ymin
corner2 = xmax, ymax
corner3 = xmin, ymax
insert!(coords, unbounded_indices[1] + 1, corner1)
insert!(coords, unbounded_indices[1] + 2, corner2)
insert!(coords, unbounded_indices[1] + 3, corner3)
elseif side1 == :right && side2 == :top
corner = xmax, ymax
insert!(coords, unbounded_indices[1] + 1, corner)
elseif side1 == :right && side2 == :left
corner1 = xmax, ymax
corner2 = xmin, ymax
insert!(coords, unbounded_indices[1] + 1, corner1)
insert!(coords, unbounded_indices[1] + 2, corner2)
elseif side1 == :right && side2 == :bottom
corner1 = xmax, ymax
corner2 = xmin, ymax
corner3 = xmin, ymin
insert!(coords, unbounded_indices[1] + 1, corner1)
insert!(coords, unbounded_indices[1] + 2, corner2)
insert!(coords, unbounded_indices[1] + 3, corner3)
elseif side1 == :top && side2 == :left
corner = xmin, ymax
insert!(coords, unbounded_indices[1] + 1, corner)
elseif side1 == :top && side2 == :bottom
corner1 = xmin, ymax
corner2 = xmin, ymin
insert!(coords, unbounded_indices[1] + 1, corner1)
insert!(coords, unbounded_indices[1] + 2, corner2)
elseif side1 == :top && side2 == :right
corner1 = xmin, ymax
corner2 = xmin, ymin
corner3 = xmax, ymin
insert!(coords, unbounded_indices[1] + 1, corner1)
insert!(coords, unbounded_indices[1] + 2, corner2)
insert!(coords, unbounded_indices[1] + 3, corner3)
elseif side1 == :left && side2 == :bottom
corner = xmin, ymin
insert!(coords, unbounded_indices[1] + 1, corner)
elseif side1 == :left && side2 == :right
corner1 = xmin, ymin
corner2 = xmax, ymin
insert!(coords, unbounded_indices[1] + 1, corner1)
insert!(coords, unbounded_indices[1] + 2, corner2)
elseif side1 == :left && side2 == :top
corner1 = xmin, ymin
corner2 = xmax, ymin
corner3 = xmax, ymax
insert!(coords, unbounded_indices[1] + 1, corner1)
insert!(coords, unbounded_indices[1] + 2, corner2)
insert!(coords, unbounded_indices[1] + 3, corner3)
end
end
end
push!(coords, coords[begin])
return coords
end

"""
get_neighbouring_boundary_edges(vor::VoronoiTessellation, e)
Expand Down Expand Up @@ -475,23 +348,28 @@ Gets the centroid of the polygon with index `i` in `vor`.
get_centroid(vor::VoronoiTessellation, i) = polygon_features(vor, i)[2]

"""
polygon_bounds(vorn::VoronoiTessellation, unbounded_extension_factor=0.0)
polygon_bounds(vorn::VoronoiTessellation, unbounded_extension_factor=0.0; include_polygon_vertices=true)
Gets the bounding box of the polygons in `vorn`. If `unbounded_extension_factor` is positive, the bounding box is extended by this factor in each direction,
proportional to the width of each axis.
If `include_polygon_vertices=true`, then the bounds both the generators and the polygons. Otherwise, only the generators
will be considered.
"""
function polygon_bounds(vorn::VoronoiTessellation, unbounded_extension_factor=0.0)
function polygon_bounds(vorn::VoronoiTessellation, unbounded_extension_factor=0.0; include_polygon_vertices=true)
F = number_type(vorn)
xmin = typemax(F)
xmax = typemin(F)
ymin = typemax(F)
ymax = typemin(F)
for i in each_polygon_vertex(vorn)
x, y = _getxy(get_polygon_point(vorn, i))
xmin = min(xmin, x)
xmax = max(xmax, x)
ymin = min(ymin, y)
ymax = max(ymax, y)
if include_polygon_vertices
for i in each_polygon_vertex(vorn)
x, y = _getxy(get_polygon_point(vorn, i))
xmin = min(xmin, x)
xmax = max(xmax, x)
ymin = min(ymin, y)
ymax = max(ymax, y)
end
end
for i in each_generator(vorn)
x, y = _getxy(get_generator(vorn, i))
Expand All @@ -500,11 +378,11 @@ function polygon_bounds(vorn::VoronoiTessellation, unbounded_extension_factor=0.
ymin = min(ymin, y)
ymax = max(ymax, y)
end
xmin -= unbounded_extension_factor * (xmax - xmin)
xmax += unbounded_extension_factor * (xmax - xmin)
ymin -= unbounded_extension_factor * (ymax - ymin)
ymax += unbounded_extension_factor * (ymax - ymin)
return xmin, xmax, ymin, ymax
_xmin = xmin - unbounded_extension_factor * (xmax - xmin)
_xmax = xmax + unbounded_extension_factor * (xmax - xmin)
_ymin = ymin - unbounded_extension_factor * (ymax - ymin)
_ymax = ymax + unbounded_extension_factor * (ymax - ymin)
return _xmin, _xmax, _ymin, _ymax
end

"""
Expand Down
89 changes: 88 additions & 1 deletion src/geometry_utils/intersections.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -197,4 +197,91 @@ function classify_and_compute_segment_intersection(a, b, c, d)
F = number_type(a)
return cert, cert_c, cert_d, (F(NaN), F(NaN))
end
end
end

"""
intersection_of_ray_with_edge(p, q, a, b)
Given a ray starting at `p` and in the direction of `q` out to
infinity, finds the intersection of the ray with the edge from `a` to `b`.
If no such intersection exists, then the coordinates of the intersection are returned
as `(NaN, NaN)`.
"""
function intersection_of_ray_with_edge(p, q, a, b)
p1, p2 = _getxy(p)
q1, q2 = _getxy(q)
a1, a2 = _getxy(a)
b1, b2 = _getxy(b)
den = a2 * q1 - a2 * p1 - b1 * p2 + b2 * p1 + b1 * q2 - b2 * q1 + a1 * (p2 - q2)
t = -(a1 * b2 - a2 * b1 - a1 * p2 + a2 * p1 + b1 * p2 - b2 * p1) / den
u = (a2 * q1 - a2 * p1 + p1 * q2 - p2 * q1 + a1 * (p2 - q2)) / den
if t < 0 || u < 0 || u > 1
return (NaN, NaN)
else
return (p1 + t * (q1 - p1), p2 + t * (q2 - p2))
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
3 changes: 2 additions & 1 deletion src/geometry_utils/polygons.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -245,4 +245,5 @@ function sort_convex_polygon!(vertices, points)
vert_to_angle = v -> (to_angle ∘ get_point)(points, v)
sort!(vertices, by=vert_to_angle)
return vertices
end
end

Loading

0 comments on commit 7ddb44a

Please sign in to comment.