[WIP] Spherical triangulation and interpolation #182
Conversation
|
|
||
| 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.
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.
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.
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.
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.
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.
I think one problem is that the mesh and convex hull here don't actually cover the full Earth - see mesh(bertin_points, faces)
There was a problem hiding this comment.
|
|
||
| faces[findall(x -> 1 in x, faces)] | ||
|
|
||
| # try NaturalNeighbours, fail miserably |
There was a problem hiding this comment.
Are there issues even when using Sibson(1) below? Or what's the misery
There was a problem hiding this comment.
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 <[email protected]>
| 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.
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.
boundary_faces is also unused
There was a problem hiding this comment.
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.
@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.
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.
On main the boundary points generate degenerate triangles, I guess I can remove those?
There was a problem hiding this comment.
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.
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.
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.
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.
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.
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.
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! |
|
The envelope would be the same, i.e. the natural neighbours from the envelope in 2D are actually the natural neighbours on the sphere. It's the actual computation of the Laplacian coordinates that could differ. So maybe this is in some way a natural neighbour interpolant. Remember that the Laplace coordinates are defined by where the sum of over the set of all natural neighbours and where |
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. |
|
That sounds right to me - spherical distance is just an inverse cosine of the dot product between Cartesian vectors, so relatively easy to compute. Circumcenters and circumradii are also easy to find, and I might end up caching them since the number of interpolation points is not necessarily very large. I don't think that this proof will extend to insertions, because we have to modify the triangulation when going back to the sphere. Specifically, since we rotate about a pivot point, we have to change the ghost vertex to point back to the pivot point. Since the ghost vertex is at infinity, inserting a point on one side of the triangulation will not change the other side, even if it should. I'm not 100% sure about this and will have to do some testing, but I suspect it is the case. I got something with the spherical Laplace approach, though there are a lot of NaNs coming through: |
|
Fixed one bug, this looks better now: |
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 |
|
Just fixed the winding order, I still have some strange looking circumcircle artifacts but this is looking a lot more reasonable than before :D |
asinghvi17
changed the title
- Caps not propagated - randsphere fixed
|
Here's a test case: 10 random points on the sphere (value 0), plus one at the north pole (value 1). 
Will probably help quite a bit with debugging... |
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