
"""
    intersect(a::Line, b::Line) -> Point

Intersection of 2 line segmens `a` and `b`.
Returns intersection_found::Bool, intersection_point
"""
function intersects(a::Line{2, T1}, b::Line{2, T2}) where {T1, T2}
    T = promote_type(T1, T2)
    v1, v2 = a; v3, v4 = b;
    MT = Mat{2,2,T,4}
    p0 = zero(Point2{T})

    verticalA = v1[1] == v2[1]
    verticalB = v3[1] == v4[1]

    # if a segment is vertical the linear algebra might have trouble
    # so we will rotate the segments such that neither is vertical
    dorotation = verticalA || verticalB

    if dorotation
        θ = T(0.0)
        if verticalA && verticalB
            θ = T(π / 4)
        elseif verticalA || verticalB # obviously true, but make it clear
            θ34 = -atan(v4[2] - v3[2], v4[1] - v3[1])
            θ12 = -atan(v2[2] - v1[2], v2[1] - v1[1])
            θ = verticalA ? θ34 : θ12
            θ = abs(θ) == T(0) ? (θ12 + θ34) / 2 : θ
            θ = abs(θ) == T(pi) ? (θ12 + θ34) / 2 : θ
        end
        rotation = MT(cos(θ), sin(θ), -sin(θ), cos(θ))
        v1 = rotation * v1
        v2 = rotation * v2
        v3 = rotation * v3
        v4 = rotation * v4
    end

    a = det(MT(v1[1] - v2[1], v1[2] - v2[2],
               v3[1] - v4[1], v3[2] - v4[2]))

    (abs(a) < eps(T)) && return false, p0 # Lines are parallel

    d1 = det(MT(v1[1], v1[2], v2[1], v2[2]))
    d2 = det(MT(v3[1], v3[2], v4[1], v4[2]))
    x = det(MT(d1, v1[1] - v2[1], d2, v3[1] - v4[1])) / a;
    y = det(MT(d1, v1[2] - v2[2], d2, v3[2] - v4[2])) / a;

    (x < prevfloat(min(v1[1], v2[1])) || x > nextfloat(max(v1[1], v2[1]))) && return false, p0
    (y < prevfloat(min(v1[2], v2[2])) || y > nextfloat(max(v1[2], v2[2]))) && return false, p0
    (x < prevfloat(min(v3[1], v4[1])) || x > nextfloat(max(v3[1], v4[1]))) && return false, p0
    (y < prevfloat(min(v3[2], v4[2])) || y > nextfloat(max(v3[2], v4[2]))) && return false, p0

    point = Point2{T}(x, y)
    # don't forget to rotate the answer back
    if dorotation
        point = transpose(rotation) * point
    end

    return true, point
end


function simple_concat(vec::AbstractVector, range, endpoint::P) where P
    result = Vector{P}(undef, length(range) + 1)
    for (i, j) in enumerate(range)
        result[i] = vec[mod1(j, length(vec))]
    end
    result[end] = endpoint
    return result
end


function consecutive_pairs(arr)
    n = length(arr)
    zip(view(arr, 1:n-1), view(arr, 2:n))
end

"""
    self_intersections(points::AbstractVector{AbstractPoint})

Finds all self intersections of polygon `points`
"""
function self_intersections(points::AbstractVector{<:AbstractPoint})
    sections = similar(points, 0)
    intersections = Int[]

    wraparound(i) = mod1(i, length(points) - 1)

    for (i, (a,b)) in enumerate(consecutive_pairs(points))
        for (j, (a2, b2)) in enumerate(consecutive_pairs(points))
            is1, is2 = wraparound(i+1), wraparound(i-1)
            if i != j && is1 != j && is2 != j && !(i in intersections) && !(j in intersections)
                intersected, p = intersects(Line(a,b), Line(a2, b2))
                if intersected
                    push!(intersections, i, j)
                    push!(sections, p)
                end
            end
        end
    end
    return intersections, sections
end

"""
    split_intersections(points::AbstractVector{AbstractPoint})

Splits polygon `points` into it's self intersecting parts. Only 1 intersection
is handled right now.
"""
function split_intersections(points::AbstractVector{<:AbstractPoint})
    intersections, sections = self_intersections(points)
    if isempty(intersections)
        return [points]
    elseif length(intersections) == 2 && length(sections) == 1
        a, b = intersections
        p = sections[1]
        a,b = min(a,b), max(a,b)
        poly1 = simple_concat(points, (a+1):(b-1), p)
        poly2 = simple_concat(points, (b+1):(length(points)+a), p)
        return [poly1, poly2]
    else
        error("More than 1 intersections can't be handled currently. Found: $intersections, $sections")
    end
end
