Intersection inconsistency between intersects and intersection for triangle and segment

[jtveiten]

Check intersection between triangle and segment, the intersection check seems to be different for intersects and intersection

julia> s1=Segment((0.5,0.5,0.), (0.5,0.5,2.0))
Segment{3,Float64}
├─ Point(0.5, 0.5, 0.0)
└─ Point(0.5, 0.5, 2.0)

julia> t1=Triangle((1.0,0.,0.),(0.,1.,0.),(0.,0.,1.))
Triangle{3,Float64}
├─ Point(1.0, 0.0, 0.0)
├─ Point(0.0, 1.0, 0.0)
└─ Point(0.0, 0.0, 1.0)

julia> intersects(s1,t1)
true

julia> intersection(s1,t1)
Intersection{Nothing}(NotIntersecting, nothing)

s1=Segment((0.5,0.5,-0.0001), (0.5,0.5,2.0))
Segment{3,Float64}
├─ Point(0.5, 0.5, -0.0001)
└─ Point(0.5, 0.5, 2.0)

julia> intersection(s1,t1)
Intersection{Point3}(Intersecting, Point(0.5, 0.5, -1.100465205072787e-17))

system info:
julia> versioninfo()
Julia Version 1.10.0
Commit 3120989f39 (2023-12-25 18:01 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: Windows (x86_64-w64-mingw32)
CPU: 8 × Intel(R) Core(TM) i5-8265U CPU @ 1.60GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, skylake)
Threads: 1 on 8 virtual cores

julia> using Pkg; Pkg.status()
Status Project.toml`
[13f3f980] CairoMakie v0.11.7
[e9467ef8] GLMakie v0.9.7
[ee78f7c6] Makie v0.20.6
[eacbb407] Meshes v0.40.4
[8913a72c] NonlinearSolve v3.5.3
[6fe1bfb0] OffsetArrays v1.13.0
[91a5bcdd] Plots v1.40.0
[1c621080] TestItems v0.1.1

[juliohm]

Thank you @jtveiten for reporting the issue. I can reproduce it.

It is a corner point intersection:

viz([s1, t1], color = ["red", "teal"])

The bug is in the intersection result. We will investigate when we find some time. Feel free to submit a PR in the meantime.

[juliohm]

@dorn-gerhard did you write our intersection(::Segment, ::Triangle) method? The fix may consist of a simple change from < to <= somewhere in the function body.

[dorn-gerhard]

Hi,
I cannot find an implementation for intersects(::Segment, ::Triangle) in src/predicates/intersects.jl
Not sure why it returns true.

No, I haven't implemented Triangle - Segment.
I was thinking about it but a while ago using a generalization of intersectparameters to more dimensions, then we got stuck in that discussion of #325 where I wanted to fix names for intersection Types. Then they got refactored which simplified a lot.

I think another refactoring might be necessary to consider intersection of objects of different dimensionality.
e.g. each object has an inner part, e.g. interval 0 < x < 1 and an outer part = boundary (segment end points, ray origin, edges of triangle, surface of 3d sphere).

• if a subset of the outer part with no inner points equals the intersection set it should be called Touching,
• if a subset of inner parts with a dimensionality less the original geometry is in the intersection set, I would call it Crossing
• if a subset of inner parts with the same dimensionality as the original geometry is in the intersection set, I would call it Overlapping (I would dismiss Pos and NegOverlapping as the result is clear from the intersection object (Segment vs. Ray)

Since we have two Geometries A and B in an intersection we could have combinations of these three cases above.
e.g.

• OverlappingTouching: A Segment that touches the edge of a Triangle
• OverlappingCrossing: A Segment that crosses a Triangle in plane yielding a segment
• TouchingCrossing: A Segment that crosses the origin of Ray (not parallel)

A Ray touching a corner of a Triangle (could be both also in the same plane) would be TouchingCrossing instead of CornerTouching
A Ray touching the edge of a Triangle (both not in the same plane) would also be called TouchingCrossing instead of EdgeTouching
A Ray overlapping with the edge of a Triangle (in the same plane) would be called OverlappingTouching
Not sure if CornerTouching and EdgeTouching is relevant for an algorithm.

Btw: segments, lines or rays that lie in the same plane as the triangle might not yield results:

t1 = Triangle((1.0,0.,0.),(0.,1.,0.),(0.,0.,1.))
s1 = Segment((1.0,0.,0.),(0.,1.,0.))
intersection(t1,s1)

[juliohm]

I cannot find an implementation for intersects(::Segment, ::Triangle) in src/predicates/intersects.jl
Not sure why it returns true.

You can use the @which macro to find the method that gets called:

julia> @which intersects(s1, t1)
intersects(c::Chain, g::Geometry)
     @ Meshes ~/.julia/dev/Meshes/src/predicates/intersects.jl:64

In this case we have a general method for Chain and Geometry. You can also discover this using the VSCode debugger by typing @enter intersects(s1, t1) and then following the steps.

I think another refactoring might be necessary.

These are very welcome. We are always doing that to improve the code organization.

Btw: segments, lines or rays that lie in the same plane as the triangle might not yield results

The intersection is the one with the bug in this issue. You can find the method here:

Meshes.jl/src/intersections/segments.jl

Line 206 in 3f1ea5e

function intersection(f, seg::Segment{3,T}, tri::Triangle{3,T}) where {T}

[dorn-gerhard]

I added an issue / suggestion: #734 with a better description for refactoring intersection types using inner and outer points and their dimensionality to define Overlapping, Crossing and Touching for general geometries

[CNOT]

I'm not sure this is the same bug, but even in the 2D case intersection/clipping return incorrect values. Here's a MWE that I hope could be helpful for rooting out the issue:

using Meshes
using CairoMakie
cg = Meshes.CartesianGrid((0.0, 0.0), (2.0, 1.0), (0.1, 0.1))
ps = Meshes.Point2[(0.0, 0.0), (1.0, 0.5), (1.0, 1.0), (2.0, 0.0)]
pol = Meshes.PolyArea(ps)
viz(pol,showfacets = true)
for i in faces(cg, 2)
    cl = i ∩ pol # Same thing happens with clip(i,pol, SutherlandHodgman())
    if cl != nothing
        viz!(i,showfacets = true, alpha = 0.5, facetcolor = :red)
    end
end
Makie.current_figure()

[juliohm]

@CNOT can you please open a separate issue with this MWE? It doesn't seem related.

[juliohm]

This bug is now fixed. The original article contains the following theorem:

There is an additional condition on the first vertex Q1 of the segment: it must not be coplanar with the triangle. In this case, we need to swap the vertices Q1 and Q2 before running the algorithm.

A patch release in on its way.