SutherlandHodgman not handling overlapping line intersections

[jnemanja]

In the clip function of the sutherlandhodgman.jl file, an intersection between lines is being pushed to a vector of points. This is fine as long as the intersection between lines is crossing and a point is returned. However, the intersection can be overlapping when a line is returned, resulting in an exception. Admittedly, this only happens after I override Meshes.atol(::Type{Float64}) = 1e-13 (which I need to avoid a different intersection problem). Still, it looks like this needs updating. I have a fix which in case of an overlapping intersection, pushes only the first point of the line. I am no quite sure if that is sufficient, I do not know if the returned line is always conveniently oriented to be able to always select the first point.

[juliohm]

Hi @jnemanja can you please share a MWE? Thanks for reporting the issue. Em sáb., 20 de jul. de 2024, 13:23, jnemanja ***@***.***> escreveu:

…

In the clip function of the sutherlandhodgman.jl file, an intersection between lines is being pushed to a vector of points. This is fine as long as the intersection between lines is crossing and a point is returned. However, the intersection can be overlapping when a line is returned, resulting in an exception. Admittedly, this only happens after I override Meshes.atol(::Type{Float64}) = 1e-13 (which I need to avoid a different intersection problem). Still, it looks like this needs updating. I have a fix which in case of an overlapping intersection, pushes only the first point of the line. I am no quite sure if that is sufficient, I do not know if the returned line is always conveniently oriented to be able to always select the first point. — Reply to this email directly, view it on GitHub <#950>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAZQW3KBYBJ4GJKETBIYSR3ZNKFJLAVCNFSM6AAAAABLGBZHUOVHI2DSMVQWIX3LMV43ASLTON2WKOZSGQZDAOJWGM3TAMA> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[jnemanja]

Here it is. Reproducible without overriding the atol.

p1 = Ngon(Point(0.0, 0.0), Point(0.0, 1.000000000000001), Point(1.0, 1.0))
p2 = Ngon(Point(0.0, 1.0), Point(0.0, 2.0), Point(1.0, 1.000000000001))
intersection(p1, p2)

[jnemanja]

I forgot to mention another part of the problem. It might even be an issue on its own, but probably it needs fixing for this issue to be properly fixed. When two overlapping lines are intersected, it just returns the first line instead of returning the line with two innermost points.

This:

l1 = Line(Point(0.0, 0.0), Point(2.0, 0.0))
l2 = Line(Point(1.0, 0.0), Point(3.0, 0.0))
get(intersection(l1, l2))

returns Line(Point(0.0, 0.0), Point(2.0, 0.0)) instead of Line(Point(1.0, 0.0), Point(2.0, 0.0)).

[juliohm]

@jnemanja regarding your example with Line, you probably need Segment:

julia> s1 = Segment((0, 0), (2, 0))
Segment
├─ Point(x: 0.0 m, y: 0.0 m)
└─ Point(x: 2.0 m, y: 0.0 m)

julia> s2 = Segment((1, 0), (3, 0))
Segment
├─ Point(x: 1.0 m, y: 0.0 m)
└─ Point(x: 3.0 m, y: 0.0 m)

julia> s1 ∩ s2
Segment
├─ Point(x: 1.0 m, y: 0.0 m)
└─ Point(x: 2.0 m, y: 0.0 m)

If you are coming from the GIS world, you will have to unlearn this incorrect naming.

I will try to take a look into your original MWE above soon.

[juliohm]

@jnemanja we can reduce your MWE to the following intersection between Lines:

julia> l1
Line
├─ a: Point(x: 0.0 m, y: 1.0 m)
└─ b: Point(x: 1.0 m, y: 1.000000000000001 m)

julia> l2
Line
├─ a: Point(x: 1.0 m, y: 1.0 m)
└─ b: Point(x: 0.0 m, y: 1.000000000000001 m)

julia> 
julia> l1 ∩ l2
Line
├─ a: Point(x: 0.0 m, y: 1.0 m)
└─ b: Point(x: 1.0 m, y: 1.000000000000001 m)

As you can see, these two lines are equal up to machine tolerance (i.e., Meshes.atol). You can find the exact algorithm that leads to this result here:

Meshes.jl/src/intersections/lines.jl

Lines 5 to 23 in 60996c4

	# The intersection type can be one of three types:
	#
	# 1. intersect at one point
	# 2. overlap at more than one point
	# 3. do not overlap nor intersect
	function intersection(f, line₁::Line, line₂::Line)
	a, b = line₁(0), line₁(1)
	c, d = line₂(0), line₂(1)
	λ₁, _, r, rₐ = intersectparameters(a, b, c, d)
	if r == rₐ == 2
	return @IT Crossing (a + λ₁ * (b - a)) f
	elseif r == rₐ == 1
	return @IT Overlapping line₁ f
	else
	return @IT NotIntersecting nothing f
	end
	end

As you can see, the algorithm is very simple to read and understand. What do you think we should do in this case?

[jnemanja]

Ah right, a silly misunderstanding on my part there for line and segment distinction. I'm not even from the GIS world, so it's a learning and not unlearning opportunity. However, does this still expose a problem in the SutherlandHodgman method? The method should use segments instead of lines, or it could produce malformed intersections.

[juliohm]

@jnemanja I commented above with more information. The Suthreland-Hodgman algorithm uses Lines. The issue you are facing is related to floating point errors.

[juliohm]

Here are a few suggestions:

1. Consider cleaning up your geometries if that is possible. We have some Repair algorithms implemented that you can check in the docs.
2. Consider testing alternative clipping algorithms such as the ones described in Add more clipping methods #578. We didn't have to implement them yet, and contributions are very welcome.

Please let me know if you need further help. I will mark the issue as invalid, given that it is not an issue of the clipping implementation. Our Zulip channel is a nice place to chat too.