Bug in Line Intersection ( GeometryBasics.intersects )

[mherreshoff]

I found a pair of lines that GeometryBasics.intersects claims don't intersect even though they do. How to reproduce:

I'm using Julia 1.7.2 with [5c1252a2] GeometryBasics v0.4.1

If you run the following code

using GeometryBasics
const GB = GeometryBasics
mkline(pair) = GB.Line(GB.Point(pair[1]...), GB.Point(pair[2]...))
bug2 =(mkline([-3.1, 15.588457268119894] => [3.1, 15.588457268119894]),
    mkline([2.0866025403784354, 17.37050807568877] => [-4.0866025403784505, 13.806406460551015]))
GB.intersects(bug2[1], bug2[2])

I get (false, [0.0, 0.0]) even though these lines should intersect. To confirm that they do I've drawn the lines using the Luxor package here: https://imgur.com/a/3sADDln

Drawing code for ease of copy paste:

using Luxor
@draw begin
    LP = Luxor.Point
    dscale=10
    setline(1)
    sethue("black")
    for (p1, p2) in bug2
        asp(p) = let 
            v = ([dscale, -dscale] .* vec(p))
            @show v
            Luxor.Point(v...)
        end
        line(asp(p1), asp(p2); action=:stroke)
    end
end

Anyway, thanks for having made the geometry package, and I'll be curious to see what the bug was.

[Drachta]

Hello,

I've encountered a similar bug.
What's weird is it's working fine in VSCode debug, but not in Julia console.

GeometryBasics v0.4.2 - Julia 1.6.5 - Windows 7

Minimal code to reproduce issue:

using GeometryBasics

line1 = Line(Point(5743.933982822018, 150.0) , Point(5885.355339059327, -50.0))

# return: intersects(line1, line2) = (false, [0.0, 0.0]) in Julia
# return: intersects(line1, line2) = (true, [5759.999999999998, 127.27922061357962]) in VSCode debug
line2 = Line(Point(5760.0, 100.0) , Point(5760.0, 140.0))
@show intersects(line1, line2)

# return: intersects(line1, line2) = (false, [0.0, 0.0]) in Julia
# return: intersects(line1, line2) = (true, [5800.0, 70.71067811865309]) in VSCode debug
line2 = Line(Point(5800.0, 140.0) , Point(5800.0, 0.0))
@show intersects(line1, line2)

[mherreshoff]

Interesting. I don't use VS Code myself, and I have to wonder whether it's loading a different version of Julia or GeometryBasics? Or maybe it's rounding the numbers somehow?

That said, it looks like both of our examples involve horizontal lines. I bet that's the edge case the GeometryBasics package is failing to deal with properly. (I was able to make the bug go away for the problem I was working on by rotating everything by one degree before checking for intersections.)

[cmcaine]

You're maybe just using different versions of GeometryBasics or Julia in vscode vs the julia console.

Turf.jl gets the right answers for each of the three examples in this issue and its implementation could be ported to this project.

julia> using Turf

julia> line(a, b, c, d) = LineString([Point([a, b]), Point([c, d])])
line (generic function with 1 method)

# Original example
julia> line_intersects(line(-3.1, 15.588457268119894, 3.1, 15.588457268119894), line(2.0866025403784354, 17.37050807568877, -4.0866025403784505, 13.806406460551015))
Point([-1.0000000000000058, 15.588457268119894])

# Drachta's examples
julia> line_intersects(line(5743.933982822018, 150.0, 5885.355339059327, -50.0), line(5760.0, 100.0, 5760.0, 140.0))
Point([5760.0, 127.27922061357884])

julia> line_intersects(line(5743.933982822018, 150.0, 5885.355339059327, -50.0), line(5800.0, 140.0, 5800.0, 0.0))
Point([5800.0, 70.71067811865477])

vs GeometryBasics#db65c41fd75f478b4079d8f4825784f307af024d:

julia> using GeometryBasics

julia> line(a, b, c, d) = Line(Point(a, b), Point(c, d))
line (generic function with 1 method)

# Original example
julia> intersects(line(-3.1, 15.588457268119894, 3.1, 15.588457268119894), line(2.0866025403784354, 17.37050807568877, -4.08660
25403784505, 13.806406460551015))
(false, [0.0, 0.0])

# Drachta's examples
# Note: x-intercept in first example is wrong here, but not in Turf.jl.
julia> intersects(line(5743.933982822018, 150.0, 5885.355339059327, -50.0), line(5760.0, 100.0, 5760.0, 140.0))
(true, [5759.999999999998, 127.27922061357958])

julia> intersects(line(5743.933982822018, 150.0, 5885.355339059327, -50.0), line(5800.0, 140.0, 5800.0, 0.0))
(true, [5800.0, 70.71067811865305])

A different issue, but I think it's nicer that Turf returns a Union{Point, Nothing} rather than a tuple with a boolean flag.

[chris-revell]

Hello, I realise this issue is quite old now, but I seem to be having this same problem: the code I have written using intersects seems to occasionally miss an intersection that is then obvious when the lines are plotted. This seems to only happen in my code when one of the lines is vertical, in contrast to the previous comments, but the MWE above also fails to find an intersection for me.
Version info is as follows (in both VSCode and REPL):

julia> versioninfo()
Julia Version 1.10.0
Commit 3120989f39b (2023-12-25 18:01 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 8 × Apple M1
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, apple-m1)
  Threads: 5 on 4 virtual cores
Environment:
  JULIA_NUM_THREADS = 4
  JULIA_EDITOR = code

pkg> status GeometryBasics
Status `~/Fibres/Project.toml`
  [5c1252a2] GeometryBasics v0.4.10

Quite perplexing!

[cmcaine]

You can use the fix I published Dec 2022: #185

Or move to a different package. My colleague and I moved to Turf.jl.

LibGEOS.jl and Turf.jl are both fine and you can use both through the GeoInterface interface, if I remember right. Meshes.jl seems to be where the development effort is now for pure-Julia geometry packages, but it's not at v1 yet.

I recommend using LibGEOS if you just want to move on to other issues because GEOS is a well-used library and probably already supports everything you want.

[chris-revell]

Ah, I did not realise that this fix had not been merged into the package. I looked at Meshes.jl but found it more fiddly than straight GeometryBasics. I'll look at LibGEOS and Turf. Thanks for the tips.

[chris-revell]

Your updated version of intersects works fine, incidentally.

[rafaqz]

We should really get that PR merged...

[chris-revell]

Is there anything we can do to help make that happen?

[cmcaine]

No, not unless the maintainer of this package asks for the PR to be modified before being merged.

[SimonDanisch]

There's not one maintainer for this repository ;)

[chris-revell]

Thanks very much!