Skip to content

Commit

Permalink
Remove ExactPredicates Dependance (#72)
Browse files Browse the repository at this point in the history 
* Update intersection point calculations

* Remove EP meet in clipping

* Fully remove meet

* Finish cleanup and fix to_points and to_edges

* Remove ExactPredicates dep

* Fix utils and comments

* Remove unused type param
  • Loading branch information
@skygering
skygering authored Mar 21, 2024
1 parent 3b6baf3 commit 4af2e75
Show file tree
Hide file tree
Showing 9 changed files with 129 additions and 176 deletions.
2 changes: 0 additions & 2 deletions Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,15 +4,13 @@ authors = ["Anshul Singhvi <[email protected]> and contributors"]
version = "0.0.1"

[deps]
ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110"
GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[compat]
ExactPredicates = "2"
GeoInterface = "1.2"
GeometryBasics = "0.4.7"
Proj = "1"
Expand Down
5 changes: 2 additions & 3 deletions src/GeometryOps.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,6 @@ using GeometryBasics
using GeometryBasics.StaticArrays
import Proj
using LinearAlgebra
import ExactPredicates
import Proj.CoordinateTransformations.StaticArrays
import Base.@kwdef

Expand All @@ -16,8 +15,8 @@ using GeoInterface.Extents: Extents
const GI = GeoInterface
const GB = GeometryBasics

const TuplePoint = Tuple{Float64,Float64}
const Edge = Tuple{TuplePoint,TuplePoint}
const TuplePoint{T} = Tuple{T, T} where T <: AbstractFloat
const Edge{T} = Tuple{TuplePoint{T},TuplePoint{T}} where T

include("primitives.jl")
include("utils.jl")
Expand Down
38 changes: 24 additions & 14 deletions src/methods/clipping/clipping_processor.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,13 +79,11 @@ function _build_a_list(::Type{T}, poly_a, poly_b) where T
b_pt1 = b_pt2
continue
end
int_pt, fracs = _intersection_point(T, (a_pt1, a_pt2), (b_pt1, b_pt2))
if !isnothing(fracs)
α, β = fracs
collinear = isnothing(int_pt)
# if no intersection point, skip this edge
if !collinear && 0 < α < 1 && 0 < β < 1
# Intersection point that isn't a vertex
# Determine if edges intersect and how they intersect
line_orient, intr1, intr2 = _intersection_point(T, (a_pt1, a_pt2), (b_pt1, b_pt2))
if line_orient != line_out # edges intersect
if line_orient == line_cross # Intersection point that isn't a vertex
int_pt, fracs = intr1
new_intr = PolyNode{T}(;
point = int_pt, inter = true, neighbor = j - 1,
crossing = true, fracs = fracs,
Expand All @@ -95,20 +93,32 @@ function _build_a_list(::Type{T}, poly_a, poly_b) where T
_add!(a_list, a_count, new_intr, n_a_edges)
push!(a_idx_list, a_count)
else
if (0 < β < 1 && (collinear || α == 0)) || (α == β == 0)
# a_pt1 is an intersection point
n_b_intrs += β == 0 ? 0 : 1
(_, (α1, β1)) = intr1
add_a1 = α1 == 0 && 0 ≤ β1 < 1
a1_β = add_a1 ? β1 : zero(T)
add_b1 = β1 == 0 && 0 < α1 < 1
b1_α = add_b1 ? α1 : zero(T)
if line_orient == line_over
(_, (α2, β2)) = intr2
if α2 == 0 && 0 ≤ β2 < 1
add_a1, a1_β = true, β2
end
if β2 == 0 && 0 < α2 < 1
add_b1, b1_α = true, α2
end
end
if add_a1
n_b_intrs += a1_β == 0 ? 0 : 1
a_list[prev_counter] = PolyNode{T}(;
point = a_pt1, inter = true, neighbor = j - 1,
fracs = fracs,
fracs = (zero(T), a1_β),
)
push!(a_idx_list, prev_counter)
end
if (0 < α < 1 && (collinear || β == 0))
# b_pt1 is an intersection point
if add_b1
new_intr = PolyNode{T}(;
point = b_pt1, inter = true, neighbor = j - 1,
fracs = fracs,
fracs = (b1_α, zero(T)),
)
a_count += 1
_add!(a_list, a_count, new_intr, n_a_edges)
Expand Down
115 changes: 83 additions & 32 deletions src/methods/clipping/intersection.jl
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,13 @@
# # Geometry Intersection
export intersection, intersection_points

"""
Enum LineOrientation
Enum for the orientation of a line with respect to a curve. A line can be
`line_cross` (crossing over the curve), `line_hinge` (crossing the endpoint of the curve),
`line_over` (colinear with the curve), or `line_out` (not interacting with the curve).
"""
@enum LineOrientation line_cross=1 line_hinge=2 line_over=3 line_out=4

"""
intersection(geom_a, geom_b, [T::Type]; target::Type)
Expand Down Expand Up @@ -118,13 +125,14 @@ function _intersection_points(::Type{T}, ::GI.AbstractTrait, a, ::GI.AbstractTra
if npoints_a > 0 && npoints_b > 0
# Loop over pairs of edges and add any intersection points to results
for i in eachindex(edges_a), j in eachindex(edges_b)
point, fracs = _intersection_point(T, edges_a[i], edges_b[j])
if !isnothing(point)
line_orient, intr1, _ = _intersection_point(T, edges_a[i], edges_b[j])
# TODO: Add in degenerate intersection points when line_over
if line_orient == line_cross || line_orient == line_hinge
#=
Determine if point is on edge (all edge endpoints excluded
except for the last edge for an open geometry)
=#
α, β = fracs
point, (α, β) = intr1
on_a_edge = (!a_closed && i == npoints_a && 0 <= α <= 1) ||
(0 <= α < 1)
on_b_edge = (!b_closed && j == npoints_b && 0 <= β <= 1) ||
Expand All @@ -138,43 +146,86 @@ function _intersection_points(::Type{T}, ::GI.AbstractTrait, a, ::GI.AbstractTra
return result
end

#= Calculates the intersection point between two lines if it exists, and as if the line
extended to infinity, and the fractional component of each line from the initial end point
to the intersection point.
Inputs:
(a1, a2)::Tuple{Tuple{::Real, ::Real}, Tuple{::Real, ::Real}} first line
(b1, b2)::Tuple{Tuple{::Real, ::Real}, Tuple{::Real, ::Real}} second line
Outputs:
(x, y)::Tuple{::Real, ::Real} intersection point
(t, u)::Tuple{::Real, ::Real} fractional length of lines to intersection
Both are ::Nothing if point doesn't exist!
Calculation derivation can be found here:
https://stackoverflow.com/questions/563198/
=#
function _intersection_point(::Type{T}, (a1, a2)::Tuple, (b1, b2)::Tuple) where T
#= Calculates the intersection points between two lines if they exists and the fractional
component of each line from the initial end point to the intersection point where α is the
fraction along (a1, a2) and β is the fraction along (b1, b2).
Note that the first return is the type of intersection (line_cross, line_hinge, line_over,
or line_out). The type of intersection determines how many intersection points there are.
If the intersection is line_out, then there are no intersection points and the two
intersections aren't valid and shouldn't be used. If the intersection is line_cross or
line_hinge then the lines meet at one point and the first intersection is valid, while the
second isn't. Finally, if the intersection is line_over, then both points are valid and they
are the two points that define the endpoints of the overlapping region between the two
lines.
Also note again that each intersection is a tuple of two tuples. The first is the
intersection point (x,y) while the second is the ratio along the initial lines (α, β) for
that point.
Calculation derivation can be found here: https://stackoverflow.com/questions/563198/ =#
function _intersection_point(::Type{T}, (a1, a2)::Edge, (b1, b2)::Edge) where T
# Return line orientation and 2 intersection points + fractions (nothing if don't exist)
line_orient = line_out
intr1 = ((zero(T), zero(T)), (zero(T), zero(T)))
intr2 = intr1
# First line runs from p to p + r
px, py = GI.x(a1), GI.y(a1)
rx, ry = GI.x(a2) - px, GI.y(a2) - py
# Second line runs from q to q + s
qx, qy = GI.x(b1), GI.y(b1)
sx, sy = GI.x(b2) - qx, GI.y(b2) - qy
# Intersection will be where p + tr = q + us where 0 < t, u < 1 and
# Intersections will be where p + αr = q + βs where 0 < α, β < 1 and
r_cross_s = rx * sy - ry * sx
Δqp_x = qx - px
Δqp_y = qy - py
point, fracs = if r_cross_s != 0 # if lines aren't parallel
t = (Δqp_x * sy - Δqp_y * sx) / r_cross_s
u = (Δqp_x * ry - Δqp_y * rx) / r_cross_s
x = px + t * rx
y = py + t * ry
(T(x), T(y)), (T(t), T(u))
elseif sx * Δqp_y == sy * Δqp_x # if parallel lines are collinear
t = (Δqp_x * rx + Δqp_y * ry) / (rx^2 + ry^2)
u = -(Δqp_x * sx + Δqp_y * sy) / (sx^2 + sy^2)
nothing, (T(t), T(u))
else
nothing, nothing
if r_cross_s != 0 # if lines aren't parallel
α = (Δqp_x * sy - Δqp_y * sx) / r_cross_s
β = (Δqp_x * ry - Δqp_y * rx) / r_cross_s
x = px + α * rx
y = py + α * ry
if 0 ≤ α ≤ 1 && 0 ≤ β ≤ 1
intr1 = (T(x), T(y)), (T(α), T(β))
line_orient = (α == 0 || α == 1 || β == 0 || β == 1) ? line_hinge : line_cross
end
elseif sx * Δqp_y == sy * Δqp_x # if parallel lines are collinear
# Determine overlap fractions
r_dot_r = (rx^2 + ry^2)
s_dot_s = (sx^2 + sy^2)
r_dot_s = rx * sx + ry * sy
b1_α = (Δqp_x * rx + Δqp_y * ry) / r_dot_r
b2_α = b1_α + r_dot_s / r_dot_r
a1_β = -(Δqp_x * sx + Δqp_y * sy) / s_dot_s
a2_β = a1_β + r_dot_s / s_dot_s
# Determine which endpoints start and end the overlapping region
n_intrs = 0
if 0 ≤ a1_β ≤ 1
n_intrs += 1
intr1 = (T.(a1), (zero(T), T(a1_β)))
end
if 0 ≤ a2_β ≤ 1
n_intrs += 1
new_intr = (T.(a2), (one(T), T(a2_β)))
n_intrs == 1 && (intr1 = new_intr)
n_intrs == 2 && (intr2 = new_intr)
end
if 0 < b1_α < 1
n_intrs += 1
new_intr = (T.(b1), (T(b1_α), zero(T)))
n_intrs == 1 && (intr1 = new_intr)
n_intrs == 2 && (intr2 = new_intr)
end
if 0 < b2_α < 1
n_intrs += 1
new_intr = (T.(b2), (T(b2_α), one(T)))
n_intrs == 1 && (intr1 = new_intr)
n_intrs == 2 && (intr2 = new_intr)
end
if n_intrs == 1
line_orient = line_hinge
elseif n_intrs == 2
line_orient = line_over
end
end
return point, fracs
return line_orient, intr1, intr2
end
91 changes: 5 additions & 86 deletions src/methods/geom_relations/geom_geom_processors.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,15 +12,6 @@ Enum for the orientation of a point with respect to a curve. A point can be
"""
@enum PointOrientation point_in=1 point_on=2 point_out=3

"""
Enum LineOrientation
Enum for the orientation of a line with respect to a curve. A line can be
`line_cross` (crossing over the curve), `line_hinge` (crossing the endpoint of the curve),
`line_over` (colinear with the curve), or `line_out` (not interacting with the curve).
"""
@enum LineOrientation line_cross=1 line_hinge=2 line_over=3 line_out=4


#=
Determines if a point meets the given checks with respect to a curve.
Expand Down Expand Up @@ -146,7 +137,7 @@ function _line_curve_process(
for j in (closed_curve ? 1 : 2):nc
c_end = _tuple_point(GI.getpoint(curve, j))
# Check if line and curve segments meet
seg_val = _segment_segment_orientation((l_start, l_end), (c_start, c_end))
seg_val, intr1, _ = _intersection_point(Float64, (l_start, l_end), (c_start, c_end))
# If segments are co-linear
if seg_val == line_over
!over_allow && return false
Expand All @@ -164,7 +155,7 @@ function _line_curve_process(
in_req_met = true
elseif seg_val == line_hinge # could cross or overlap
# Determine location of intersection point on each segment
(α, β) = _find_intersect_fracs(l_start, l_end, c_start, c_end)
(_, (α, β)) = intr1
if ( # Don't consider edges of curves as they can't cross
(!closed_line && ((α == 0 && i == 2) || (α == 1 && i == nl))) ||
(!closed_curve && ((β == 0 && j == 2) || (β == 1 && j == nc)))
Expand All @@ -180,7 +171,8 @@ function _line_curve_process(
α, β, l_start, l_end, c_start, c_end,
i, line, j, curve,
)
if _segment_segment_orientation(l, c) == line_hinge
next_val, _, _ = _intersection_point(Float64, l, c)
if next_val == line_hinge
!cross_allow && return false
else
!over_allow && return false
Expand Down Expand Up @@ -221,30 +213,6 @@ function _find_new_seg(i, ls, le, cs, ce)
return i, ls, break_off
end

#= Find where line and curve segments intersect by fraction of length. α is the fraction of
the line (ls to le) and β is the traction of the curve (cs to ce). All inputs are tuples. =#
function _find_intersect_fracs(ls, le, cs, ce)
point, fracs = _intersection_point(
Float64,
(ls, le),
(cs, ce)
)
(α, β) = if !isnothing(point)
fracs
else # line and curve segments are parallel
if equals(ls, cs)
(0, 0)
elseif equals(ls, ce)
(0, 1)
elseif equals(le, cs)
(1, 0)
else # equals(l_end, c_end)
(1, 1)
end
end
return α, β
end

#= Find next set of segments needed to determine if given hinge segments cross or not.=#
function _find_hinge_next_segments(α, β, ls, le, cs, ce, i, line, j, curve)
next_seg = if β == 1
Expand Down Expand Up @@ -537,52 +505,6 @@ function _point_filled_curve_orientation(
return iseven(k) ? out : in
end

#=
Determines the type of interaction between two line segments. If the segments
`cross`, this means that they have a single intersection point that isn't on
either of their enpoints. If they form a `hinge`, they meet at one of the
segments endpoints. If they are `over`, then they are co-linear for at least
some of the length of the segments. Finally, if they are `out`, then the
segments are disjoint.
Point should be an object of point trait and curve should be an object with a
linestring or linearring trait.
Can provide values of in, on, and out keywords, which determines return values
for each scenario.
=#
function _segment_segment_orientation(
(a_point, b_point), (c_point, d_point);
cross::T = line_cross, hinge::T = line_hinge,
over::T = line_over, out::T = line_out,
) where T
(ax, ay) = Float64.(a_point)
(bx, by) = Float64.(b_point)
(cx, cy) = Float64.(c_point)
(dx, dy) = Float64.(d_point)
meet_type = ExactPredicates.meet((ax, ay), (bx, by), (cx, cy), (dx, dy))
# Lines meet at one point within open segments
meet_type == 1 && return cross
# Lines don't meet at any points
meet_type == -1 && return out
# Lines meet at one or more points within closed segments
if _isparallel(((ax, ay), (bx, by)), ((cx, cy), (dx, dy)))
min_x, max_x = cx < dx ? (cx, dx) : (dx, cx)
min_y, max_y = cy < dy ? (cy, dy) : (dy, cy)
if (
((ax ≤ min_x && bx ≤ min_x) || (ax ≥ max_x && bx ≥ max_x)) &&
((ay ≤ min_y && by ≤ min_y) || (ay ≥ max_y && by ≥ max_y))
)
# a_point and b_point are on the same side of segment, don't overlap
return hinge
else
return over
end
end
# if lines aren't parallel then they must hinge
return hinge
end

#=
Determines the types of interactions of a line with a filled-in curve. By
filled-in curve, I am referring to the exterior ring of a poylgon, for example.
Expand Down Expand Up @@ -636,10 +558,7 @@ function _line_filled_curve_interactions(
for j in 1:nc
c_end = _tuple_point(GI.getpoint(curve, j))
# Check if two line and curve segments meet
seg_val = _segment_segment_orientation(
(l_start, l_end),
(c_start, c_end),
)
seg_val, _, _ = _intersection_point(Float64, (l_start, l_end), (c_start, c_end))
if seg_val != line_out
# If line and curve meet, then at least one point is on boundary
on_curve = true
Expand Down
Loading

0 comments on commit 4af2e75

Please sign in to comment.