Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Sg/clipping errors #146

Merged
merged 43 commits into from
May 29, 2024
Merged

Sg/clipping errors #146

merged 43 commits into from
May 29, 2024

Conversation

skygering
Copy link
Collaborator

In this PR, we re-add ExactPredicates through the Predicates file. In this file, I have made "exact" and "non-exact" versions of each predicate, which correspond to the use to ExactPredicates for the calculations vs not using ExactPredicates. I have added an exact keyword to denote which predicate should be used throughout the rest of the functions that depend on these predicates. These exact keywords don't make it up to the user level, but rather stay one level below that within helper functions.

All clipping uses exact = true and the geometry relations use exact = false. This is because something in the geometry relation code causes some of the tests to fail when exact = true. I think that has something to do with line 621 in geom_geom_processors.jl which checks the midpoint of a line to see if it is in, out, or on a filled curve. I think then when it is supposed to be "on", due to floating point issues when adding two points together and dividing by two to get the midpoint, it is reading as "in" or "out" (however, this is another issue / PR). On the other hand, clipping needs exact = true for quite a few edge cases.

In my opinion, I don't think we should keep the exact keyword in the long term. ExactPredicates should be improved to use adaptive predicates as well and then all of the code should adjust itself to be as exact as needed to get the most correct answer.

However, for now, the exact keyword allowed me to update the clipping without breaking the current functionality for the geometry relations. It might also be a while before we have the adaptive solutions working, so it is nice to have this middle ground.

I also re-wrote the intersection_point function to make it much more intuitive.

@asinghvi17 also added some benchmarking plots.

Another thing to open an issue about is checking how coverage and cuts deal with some of these edge cases. I haven't really thought too hard about how they will interact with some of the changes I have made. Also, the intersection_points function. It isn't used anywhere right now, but it doesn't currently output the intersection points from collinear segments. So this needs to just be re-done. I think this will eventually be used within clipping if we want to optimize the finding of intersection points. So maybe that should be a combined PR.

skygering and others added 30 commits April 20, 2024 10:21
Just so we have the code somewhere + can run it after changes
so that we can override later if necessary.  Plus, this preserves hygiene.
Add benchmark + refer to Predicates directly
@skygering
Copy link
Collaborator Author

skygering commented May 22, 2024

Some examples from my work using GeometryOps!

At all timesteps, the colliding polygons' (ice floes!) intersections are calculated. This simulation is 10,000 timesteps and starts with 500 polygons.

shear_flow.mp4

Every 150 timesteps, the polygons are welded together, calculating unions. This simulation is 5,000 timesteps and starts with 500 polygons. It is driven by compression walls that move inward to "squish" the polygons together and are themselves polygons. This isn't a very physical simulation and rather to push the bounds on the geometry. They bounce around because they don't actually fit in the remaining space so the forces from the invisible walls become pretty big.

packed_welding.mp4

Every 150 timesteps the floes "ridge" and "raft", calculating differences, between both the ice polygons and the wall polygons. This simulation has the same setup as the welding one above.

ridging_rafting.mp4

curve,
)
p_end = i ≤ npoints ? _tuple_point(ipoints[i]) : l_end
mid_val = _point_filled_curve_orientation((p_start .+ p_end) ./ 2, curve; exact)
Copy link
Member

@asinghvi17 asinghvi17 May 22, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You would have to perform that division and addition within the predicate, and add a grouping statement so that the error checking can handle it!

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think I can actually get rid of the checks when it is checking the midpoint of a shared edge by re-writing the logic a bit! We still might need the increased precision for other cases, but that is the one I see the greatest issues with.

@skygering
Copy link
Collaborator Author

skygering commented May 23, 2024

I debugged another clipping bug and also two simplification bugs that popped up when running my large simulations.

How do we feel about the PR? Good to merge? Also, I feel like we should do a patch release since these are all basically bugs. If so, I will open the issues mentioned above in the PR description.

@rafaqz
Copy link
Member

rafaqz commented May 24, 2024

Sorry I haven't had time to review, pretty busy currently (but happy for you to merge)

@skygering skygering requested a review from asinghvi17 May 28, 2024 15:42
Copy link
Member

@asinghvi17 asinghvi17 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good to me! Had a couple of minor comments. Documentation can be added in the next release if this is required for local work.

Comment on lines 45 to 51
c_val = if c_t1 ≈ c_t2
0
else
c = c_t1 - c_t2
c > 0 ? 1 : -1
end
return c_val
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
c_val = if c_t1 c_t2
0
else
c = c_t1 - c_t2
c > 0 ? 1 : -1
end
return c_val
return sign(c_t1 - c_t2)

maybe? It does the same thing except that sign(-0.0) == -0.0 (but -0.0 == 0).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hm won't this change the isapprox behavior?

Copy link
Member

@asinghvi17 asinghvi17 May 29, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How do you mean? Arithmetically they are the same, and -0.0 is not less than 0...

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am checking if isappox(ct1, ct2) first, but I think this could replace code in the else statement.

src/methods/geom_relations/coveredby.jl Show resolved Hide resolved
src/methods/geom_relations/disjoint.jl Outdated Show resolved Hide resolved
@skygering skygering merged commit 1474f1d into JuliaGeo:main May 29, 2024
3 checks passed
@skygering skygering deleted the sg/clipping_errors branch May 29, 2024 21:27
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants