-
Notifications
You must be signed in to change notification settings - Fork 4
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
[WIP] Spherical triangulation and interpolation #182
base: main
Are you sure you want to change the base?
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Just some comments about API. Nothing too important since this is just an experiment anyway
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Just a quick look for now. Did you have any luck trying to use ghost vertices?
|
||
ch = DelTri.convex_hull(projected_points) # assumes each point is in the triangulation | ||
boundary_nodes = DelTri.get_vertices(ch) | ||
bertin_boundary_poly = GI.Polygon([GI.LineString(DelTri.get_points(ch)[DelTri.get_vertices(ch)])]) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What do you need bertin_boundary_poly
for? Are you just using it to see if a point is inside the triangulation? You could look at DelaunayTriangulation.dist(tri, q)
I think
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I was trying to see if I could filter out points that are outside the convex hull and therefore not within any triangle of tri
, but no luck with that.
mat = [ | ||
if GO.contains(bertin_boundary_poly, (x, y)) | ||
itp(x, y; method = Nearest()) | ||
else | ||
NaN | ||
end | ||
for (x, y) in bertin_points | ||
] | ||
# TODO: this currently doesn't work, because some points are not inside a triangle and so cannot be located. | ||
# Options are: | ||
# 1. Reject all points outside the convex hull of the projected points. (Tried but failed) | ||
# 2. ... |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What goes wrong here? Can you just do something like
xv, yv = first.(bertin_points), last.(bertin_points)
vals = itp(xv, yv, method = Nearest(), project = false) # Nearest() doesn't seem to use project...
vals[identify_exterior_points(xv, yv, itp)] .= NaN
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
ERROR: TaskFailedException
nested task error: The point, (-9.986746533293853e6, -6.557879586899938e6), could not be located.
Stacktrace:
[1] brute_force_search(tri::DelaunayTriangulation.Triangulation{…}, q::Tuple{…}; itr::Set{…})
@ DelaunayTriangulation ~/.julia/dev/DelaunayTriangulation/src/algorithms/point_location/brute_force.jl:36
[2] brute_force_search
@ ~/.julia/dev/DelaunayTriangulation/src/algorithms/point_location/brute_force.jl:31 [inlined]
[3] restart_jump_and_march(tri::DelaunayTriangulation.Triangulation{…}, q::Tuple{…}, store_history::Val{…}, history::Nothing, rng::Random.TaskLocalRNG, maxiters::Int64, cur_iter::Int64, concavity_protection::Bool, num_restarts::Int64, use_barriers::Val{…})
@ DelaunayTriangulation ~/.julia/dev/DelaunayTriangulation/src/algorithms/point_location/jump_and_march.jl:1236
[4] _jump_and_march(tri::DelaunayTriangulation.Triangulation{…}, q::Tuple{…}, k::Int64, store_history::Val{…}, history::Nothing, rng::Random.TaskLocalRNG, maxiters::Int64, cur_iter::Int64, concavity_protection::Bool, num_restarts::Int64, use_barriers::Val{…})
@ DelaunayTriangulation ~/.julia/dev/DelaunayTriangulation/src/algorithms/point_location/jump_and_march.jl:1145
[5] restart_jump_and_march(tri::DelaunayTriangulation.Triangulation{…}, q::Tuple{…}, store_history::Val{…}, history::Nothing, rng::Random.TaskLocalRNG, maxiters::Int64, cur_iter::Int64, concavity_protection::Bool, num_restarts::Int64, use_barriers::Val{…})
@ DelaunayTriangulation ~/.julia/dev/DelaunayTriangulation/src/algorithms/point_location/jump_and_march.jl:1234
--- the last 2 lines are repeated 23 more times ---
[52] _jump_and_march(tri::DelaunayTriangulation.Triangulation{…}, q::Tuple{…}, k::Int64, store_history::Val{…}, history::Nothing, rng::Random.TaskLocalRNG, maxiters::Int64, cur_iter::Int64, concavity_protection::Bool, num_restarts::Int64, use_barriers::Val{…})
@ DelaunayTriangulation ~/.julia/dev/DelaunayTriangulation/src/algorithms/point_location/jump_and_march.jl:1145
[53] #jump_and_march#87
@ ~/.julia/dev/DelaunayTriangulation/src/algorithms/point_location/jump_and_march.jl:739 [inlined]
[54] jump_to_voronoi_polygon(tri::DelaunayTriangulation.Triangulation{…}, q::Tuple{…}; kwargs::@Kwargs{…})
@ DelaunayTriangulation ~/.julia/dev/DelaunayTriangulation/src/algorithms/point_location/nearest_neighbour.jl:24
[55] jump_to_voronoi_polygon
@ ~/.julia/dev/DelaunayTriangulation/src/algorithms/point_location/nearest_neighbour.jl:23 [inlined]
[56] _compute_nearest_coordinates(tri::DelaunayTriangulation.Triangulation{…}, interpolation_point::Tuple{…}, cache::NaturalNeighbours.NaturalNeighboursCache{…}; kwargs::@Kwargs{…})
@ NaturalNeighbours ~/.julia/packages/NaturalNeighbours/HJ2P0/src/interpolation/coordinates/nearest.jl:10
[57] _compute_nearest_coordinates
@ ~/.julia/packages/NaturalNeighbours/HJ2P0/src/interpolation/coordinates/nearest.jl:1 [inlined]
[58] #compute_natural_coordinates#28
@ ~/.julia/packages/NaturalNeighbours/HJ2P0/src/interpolation/coordinates/nearest.jl:19 [inlined]
[59] compute_natural_coordinates
@ ~/.julia/packages/NaturalNeighbours/HJ2P0/src/interpolation/coordinates/nearest.jl:18 [inlined]
[60] #_eval_interp#21
@ ~/.julia/packages/NaturalNeighbours/HJ2P0/src/interpolation/eval.jl:42 [inlined]
[61] _eval_interp
@ ~/.julia/packages/NaturalNeighbours/HJ2P0/src/interpolation/eval.jl:40 [inlined]
[62] (::NaturalNeighbours.NaturalNeighboursInterpolant{…})(x::Float64, y::Float64, id::Int64; parallel::Bool, method::Nearest{…}, kwargs::@Kwargs{…})
@ NaturalNeighbours ~/.julia/packages/NaturalNeighbours/HJ2P0/src/interpolation/interpolate.jl:187
[63] NaturalNeighboursInterpolant
@ ~/.julia/packages/NaturalNeighbours/HJ2P0/src/interpolation/interpolate.jl:180 [inlined]
[64] macro expansion
@ ~/.julia/packages/NaturalNeighbours/HJ2P0/src/interpolation/interpolate.jl:204 [inlined]
[65] (::NaturalNeighbours.var"#71#threadsfor_fun#18"{…})(tid::Int64; onethread::Bool)
@ NaturalNeighbours ./threadingconstructs.jl:215
[66] #71#threadsfor_fun
@ ./threadingconstructs.jl:182 [inlined]
...and 7 more exceptions.
Stacktrace:
[1] threading_run(fun::NaturalNeighbours.var"#71#threadsfor_fun#18"{…}, static::Bool)
@ Base.Threads ./threadingconstructs.jl:172
[2] macro expansion
@ ./threadingconstructs.jl:220 [inlined]
[3] (::NaturalNeighbours.NaturalNeighboursInterpolant{…})(vals::Vector{…}, x::Vector{…}, y::Vector{…}; parallel::Bool, method::Nearest{…}, kwargs::@Kwargs{…})
@ NaturalNeighbours ~/.julia/packages/NaturalNeighbours/HJ2P0/src/interpolation/interpolate.jl:202
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
That's... concerning. That should be impossible since it should be finding it in a ghost triangle at the very least... will take a look and actually run the code later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
|
||
faces[findall(x -> 1 in x, faces)] | ||
|
||
# try NaturalNeighbours, fail miserably |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Are there issues even when using Sibson(1)
below? Or what's the misery
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yeah I'm not quite sure what the issue is here, in the worst case I can always sample in stereographic space but that seems pretty bad...
Co-authored-by: Daniel VandenHeuvel <95613936+DanielVandH@users.noreply.github.com>
boundary_points = reverse([ | ||
(-far_value, -far_value), | ||
(-far_value, far_value), | ||
(far_value, far_value), | ||
(far_value, -far_value), | ||
(-far_value, -far_value), | ||
]) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is unused. Were there issues? You mentioned an issue with convert_boundary_points_to_indices
, is it related?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
boundary_faces
is also unused
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The issue here was that a lot of faces were connecting to two boundary points and one data point, which is effectively a degenerate face and thus useless.
|
||
projected_points = GO.reproject(GO.tuples(points), source_crs = GFT.EPSG(4326), target_crs = "+proj=bertin1953 +lon_0=-16.5 +lat_0=-42 +x_0=7.93 +y_0=0.09") | ||
|
||
ch = DelTri.convex_hull(projected_points) # assumes each point is in the triangulation |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@asinghvi17 projected_points
is three-dimensional. The convex_hull
function is treating it by ignoring the z
coordinate of each point. What is supposed to be the input here? This is the issue you're having with NaturalNeighbours - the triangulation is very invalid.
julia> DelTri.validate_triangulation(tri) # on main
The adjacent2vertex map is inconsistent with the adjacent map. The vertex 1953 and edge (2026, 2008) are inconsistent with the associated edge set Set([(2026, 2008), (1951, 2026), (2008, 1952), (1952, 1951)]) from the adjacent2vertex map, as (2026, 2008) maps to -1 instead of 1953.
The vertex 2365, which defines the edge (-1, 2365) with adjacent vertex 2365, is inconsistent with the graph which fails to include at least one of 2037 and -1 in 2365's neighbourhood.
The triangle (2379, 2381, 2382) is incorrectly stored in the adjacent map, with the edge (2382, 2379) instead mapping to -1.
The Delaunay criterion does not hold for the triangle-vertex pair ((1879, 763, 1464), 1444).
The graph is inconsistent with the triangle set. The triangle (1732, 2039, -1) is in the triangle set but either 2039 or -1 are not in 1732's neighbourhood.
The triangle (757, 691, 136) is negatively oriented.
The edge iterator has duplicate edges.
The solid_edge iterator has duplicate edges.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
projected_points
' z coordinate only carries metadata, so it's fine for it to be three-dimensional. That being said I haven't checked the triangulation in stereographic space, it might have issues there too.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
On main the boundary points generate degenerate triangles, I guess I can remove those?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm surprised that it's generating triangles where some of the vertices repeat..
The issue to me just looks like the triangles are somehow calculated incorrectly. In the plot I show at least the triangulation is non-planar
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
That's not entirely unexpected, since the triangulation is computed in stereographic space (where the triangulation generated is planar, correct, etc) and not in the projected space that we are currently in.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, that would be a problem. Maybe it's best to create some kind of spherical triangulation type, which can also handle spherical Voronoi? That sounds like a lot of effort though.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I don't really know how you would get around this issue. It just seems like a bug somewhere to me.. the projection isn't an isometry but it should be conformal, right? I think conformal maps would show up better than what we see in the image above, e.g. it should preserve angles. The crossings seem to come only from boundary vertices (from a visual check anyway, who knows about the interior) which suggests you're missing something when handling only the interior triangles.
If everything is fine and there are in fact no bugs in the actual creation and my thinking is incorrect, it would be interesting to know how the other library handles this. Maybe there are things I need to do on DelaunayTriangulation's side to help
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ah, to clarify - the triangulation is done in stereographic space, where the triangulation would be completely valid. I'm trying to interpolate in the Bertin 1953 projection, which doesn't have the same conformal property as the stereographic projection. I guess what I was trying to do is actually the incorrect thing, and I should be figuring out how to interpolate on the sphere instead...
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I see.. my thinking was that you can do everything on the valid triangulation and then lift it back up to the sphere. I'll probably have to look more at the algorithm itself to understand.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
At least for interpolation, I would expect that to have issues when getting the values of points near the chosen pivot point in spherical space - there's always a point which wouldn't be connected to the triangulation in stereographic space, and neighbours in spherical space would be separated by almost the entire domain of points.
The first problem could be solved by injecting a rectangular boundary, but I'm not really sure how to deal with the second problem.
@DanielVandH is it possible to re-implement Laplace coordinates if I simply insert points into the triangulation in stereographic space? Consider the following workflow:
|
That would probably work. Might not be truly natural neighbour interpolation from the sphere (I think? I don't know its definition on the sphere exactly) but it would be good enough I reckon. Let me know if it's not clear how to use some of the functions in NaturalNeighbours.jl for this. I can also make |
Yeah, that's what I was afraid of :(. There might be a way to mimic the part of the Bowyer-Watson algorithm we need on the sphere, so I might be able to reimplement natural neighbour implementation on the sphere with a few minor adjustments. Will give implementing that a try and report back shortly! |
Just as an actual proof of this in case you don't know the full details: Consider adding a point into the triangulation on a sphere. To insert this point, you would project the sphere, insert it into the 2D triangulation, and then go back to the sphere. The natural neighbours are the vertices that share an edge with this new point. Since the projection doesn't change the connectivities, these are the same vertices as in 2D. i.e. the natural neighbours are exactly the same. |
Ah. You're probably right.. wasn't thinking about the infinity issues. It'd be interesting to see how natural neighbour interpolation on a sphere is actually done (I don't have access to that book you asked for on Slack either) to know what to do. I do some NaN handling in some places in NaturalNeighbours.jl I think. Might want to search for |
- Caps not propagated - randsphere fixed
This works, but the code is a bit messy and I implemented a bunch of things myself. Would love to return a custom point type for spherical coordinates that encodes Cartesian coords on the unit sphere, but that would be a bit much for an experiment. Maybe when we actually do it :D
Screen.Recording.2024-07-21.at.10.52.39.PM_1.mp4