diff --git a/Project.toml b/Project.toml index 5c8df722..4c7183ff 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GeoMakie" uuid = "db073c08-6b98-4ee5-b6a4-5efafb3259c6" authors = ["Makie.jl Contributors"] -version = "0.7.5" +version = "0.7.6" [deps] Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" diff --git a/src/GeoMakie.jl b/src/GeoMakie.jl index ccba8c6b..7d82806f 100644 --- a/src/GeoMakie.jl +++ b/src/GeoMakie.jl @@ -45,6 +45,7 @@ const Text = Makie.Text Base.convert(::Type{Rect{N, Float64}}, x::Rect{N}) where N = Rect{N, Float64}(x) include("makie_piracy.jl") +include("triangulation3d.jl") include("geojson.jl") # GeoJSON/GeoInterface support include("conversions.jl") include("data.jl") diff --git a/src/triangulation3d.jl b/src/triangulation3d.jl new file mode 100644 index 00000000..42754721 --- /dev/null +++ b/src/triangulation3d.jl @@ -0,0 +1,113 @@ +using Makie.LinearAlgebra +using GeometryBasics +using GeometryOps.ExactPredicates +using GeometryOps.ExactPredicates.Codegen +using GeometryOps.ExactPredicates.StaticArrays: SVector + +# This function is an "exact predicate" that is robust to floating point error +# May not be necessary, but it's here if we need it. +@genpredicate function _collinear(p :: 3, q :: 3, r :: 3) + pq = q - p + pr = r - p + Codegen.group!(pq...) + Codegen.group!(pr...) + # Cross product of vectors will be zero if points are collinear + cross = SVector{3}( + pq[2]*pr[3] - pq[3]*pr[2], + pq[3]*pr[1] - pq[1]*pr[3], + pq[1]*pr[2] - pq[2]*pr[1] + ) + return ExactPredicates.inp(cross, cross) # Will be 0 if collinear, positive otherwise +end + +# This should be removed once https://github.com/MakieOrg/Makie.jl/pull/4584 is merged! +# TODO + +# We don't want to override the general case poly_convert for all polygons, because that's piracy, +# but we _can_ override it for the specific case of a 3D polygon that is being transformed by a function +# that is a subtype of Union{<: Proj.Transformation, <: GeoMakie.Geodesy.ECEFfromLLA} +function Makie.poly_convert(polygon::GeometryBasics.Polygon, transform_func::Union{<: Proj.Transformation, <: GeoMakie.Geodesy.ECEFfromLLA}) + + outer = GeometryBasics.metafree(GeometryBasics.coordinates(polygon.exterior)) + PT = Makie.float_type(outer) + points = [Makie.apply_transform(transform_func, outer)] + points_flat = PT[outer;] + for inner in polygon.interiors + inner_points = GeometryBasics.metafree(GeometryBasics.coordinates(inner)) + append!(points_flat, inner_points) + push!(points, Makie.apply_transform(transform_func, inner_points)) + end + + # Shortcut if the transformation is 2D -> 2D + if points isa Vector{Vector{<: Makie.VecTypes{2}}} + faces = GeometryBasics.earcut_triangulate(points) + return GeometryBasics.Mesh(points_flat, faces) + end + + # We assume the polygon lies on a plane, and thus seek to find that plane, + # so we can use it to project the polygon into 2D, and then call earcut_triangulate + # on the projected polygon. + + # First, we extract three unique and independent (non-collinear) points from the polygon. + p1, p2, p3 = extract_three_unique_and_independent_points(points) + + # Now, we can find a plane from these points: + + # Define a plane that can be used to project the polygon into 2D + v1 = p2 - p1 + v2 = p3 - p1 + normal = cross(v1, v2) + + # `x` and `y` are the vectors that define the plane. + x = v1 + y = cross(normal, x) + + + # Project the polygon into 2D + projected_polygon = map(ring -> map(p -> Point2{Float64}(dot(p, x), dot(p, y)), ring), points) + + # Now, call earcut_triangulate on the projected polygon, which is 2D + faces = GeometryBasics.earcut_triangulate(projected_polygon) + return GeometryBasics.Mesh(points_flat, faces) +end + +function extract_three_unique_and_independent_points(points::Vector{Vector{PT}}) where PT <: Makie.VecTypes + p1, p2, p3 = points[1][1], points[1][2], points[1][3] + + if p1 == p2 || p1 == p3 || p2 == p3 + if length(points[1]) <= 3 + error("Polygon has only three points and they are all the same, we can't triangulate this!") + elseif p1 == p2 == p3 + new_point_idx = findfirst(p -> p != p1, points[1]) + if isnothing(new_point_idx) + error("All points in the polygon are the same, we can't triangulate this!") + end + p2 = points[1][new_point_idx] + new_point_idx = findfirst(p -> p != p1 && p != p2, points[1]) + if isnothing(new_point_idx) + error("Only found two unique points in the polygon, we can't triangulate this!") + end + p3 = points[1][new_point_idx] + elseif p1 == p2 + p2 = points[1][4] + elseif p1 == p3 + p3 = points[1][4] + elseif p2 == p3 + p2 = points[1][4] + end + end + + # Account for collinear points + if _collinear(Point3{Float64}(p1), Point3{Float64}(p2), Point3{Float64}(p3)) == 0 # collinear, all the points lie on the same line + if length(points[1]) <= 3 + error("Polygon has only three points and they are all collinear, we can't triangulate this!") + end + new_point_idx = findfirst(p -> _collinear(Point3{Float64}(p1), Point3{Float64}(p2), Point3{Float64}(p)) != 0, points[1]) + if isnothing(new_point_idx) + error("All points in the polygon are collinear, we can't triangulate this!") + end + p3 = points[1][new_point_idx] + end + + return p1, p2, p3 +end diff --git a/test/triangulation3d.jl b/test/triangulation3d.jl new file mode 100644 index 00000000..3697875e --- /dev/null +++ b/test/triangulation3d.jl @@ -0,0 +1,106 @@ +using Test +using GeoMakie, Makie +using GeometryBasics +import GeometryOps as GO, GeoInterface as GI +using Geodesy +using LinearAlgebra + +@testset "3D Polygon Triangulation" begin + @testset "Simple planar polygon" begin + # Create a simple square in 3D space + points = [ + Point3f(0, 0, 0), + Point3f(1, 0, 0), + Point3f(1, 1, 0), + Point3f(0, 1, 0) + ] + poly = Polygon(points) + + # Test triangulation + mesh = Makie.poly_convert(poly) + faces = GeometryBasics.faces(mesh) + @test length(faces) == 2 # A square should decompose into 2 triangles + @test faces isa Vector{GeometryBasics.GLTriangleFace} + end + + @testset "Rotated planar polygon" begin + # Create a square rotated 45° around y-axis + points = [ + Point3f(0, 0, 0), + Point3f(1/√2, 0, 1/√2), + Point3f(1/√2, 1, 1/√2), + Point3f(0, 1, 0) + ] + poly = Polygon(points) + + mesh = Makie.poly_convert(poly) + faces = GeometryBasics.faces(mesh) + @test length(faces) == 2 + end + + @testset "Error cases" begin + # Test collinear points + collinear_points = [ + Point3f(0, 0, 0), + Point3f(1, 1, 1), + Point3f(2, 2, 2) + ] + @test_throws ErrorException Makie.poly_convert(Polygon(collinear_points)) + + # Test duplicate points + duplicate_points = [ + Point3f(0, 0, 0), + Point3f(0, 0, 0), + Point3f(0, 0, 0) + ] + @test_throws ErrorException Makie.poly_convert(Polygon(duplicate_points)) + end + + @testset "Complex spherical polygon" begin + # Create a test case similar to your diagnostic code + londs = [0.0, 90.0, 180.0, 270.0] + latds = [45.0, 45.0, 45.0, 45.0] + + # Convert to 3D points using your transformation + transf = GeoMakie.Geodesy.ECEFfromLLA(GeoMakie.Geodesy.WGS84()) + points = [Point2(λ, φ) for (λ, φ) in zip(londs, latds)] + poly = Polygon(points) + + # Transform to 3D + transformed_poly = GO.transform(poly) do point + Makie.apply_transform(transf, point) + end |> x -> GO.GI.convert(GeometryBasics, x) + + mesh = Makie.poly_convert(transformed_poly) + faces = GeometryBasics.faces(mesh) + + meshfrom2d = Makie.poly_convert(poly, transf) + facesfrom2d = GeometryBasics.faces(meshfrom2d) + + @test faces == facesfrom2d + + @test length(faces) > 0 # Should produce at least one triangle + @test faces isa Vector{GeometryBasics.GLTriangleFace} + end + + @testset "Polygon with interior" begin + # Create a square with a square hole + outer = [ + Point3f(0, 0, 0), + Point3f(2, 0, 0), + Point3f(2, 2, 0), + Point3f(0, 2, 0) + ] + inner = [ + Point3f(0.5, 0.5, 0), + Point3f(1.5, 0.5, 0), + Point3f(1.5, 1.5, 0), + Point3f(0.5, 1.5, 0) + ] + + poly = Polygon(outer, [inner]) + mesh = Makie.poly_convert(poly) + faces = GeometryBasics.faces(mesh) + @test length(faces) > 4 # Should need more than 4 triangles to fill the ring + end +end \ No newline at end of file