Add MaxLengthDiscretization (#1068)

* Add 'MaxLengthDiscretization'

* Add more methods

* Update src/discretization/maxlength.jl

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

* Update code

* Apply suggestions

* Update code

* Format

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

* Add tests

---------

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

Author: eliascarv

src/Meshes.jl

@@ -506,6 +506,7 @@ export
   DelaunayTriangulation,
   ManualDiscretization,
   RegularDiscretization,
+  MaxLengthDiscretization,
   discretize,
   discretizewithin,
   simplexify,

src/discretization.jl

@@ -190,8 +190,6 @@ simplexify(box::Box{𝔼{2}}) = discretize(box, FanTriangulation())
 
 simplexify(box::Box{𝔼{3}}) = discretize(box, ManualDiscretization())
 
-simplexify(seg::Segment) = SimpleMesh(pointify(seg), GridTopology(1))
-
 function simplexify(chain::Chain)
   np = nvertices(chain) + isclosed(chain)
   ip = isperiodic(chain)
@@ -263,3 +261,4 @@ include("discretization/held.jl")
 include("discretization/delaunay.jl")
 include("discretization/manual.jl")
 include("discretization/regular.jl")
+include("discretization/maxlength.jl")

src/discretization/maxlength.jl

@@ -0,0 +1,58 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    MaxLengthDiscretization(length)
+
+Discretize geometries into parts with sides of maximum `length` in length units (default to meters).
+"""
+struct MaxLengthDiscretization{ℒ<:Len} <: DiscretizationMethod
+  length::ℒ
+  MaxLengthDiscretization(length::ℒ) where {ℒ<:Len} = new{float(ℒ)}(length)
+end
+
+MaxLengthDiscretization(length) = MaxLengthDiscretization(addunit(length, u"m"))
+
+function discretize(box::Box, method::MaxLengthDiscretization)
+  sizes = ceil.(Int, _sides(box) ./ method.length)
+  discretize(box, RegularDiscretization(sizes))
+end
+
+function discretize(segment::Segment, method::MaxLengthDiscretization)
+  size = ceil(Int, _measure(segment) / method.length)
+  discretize(segment, RegularDiscretization(size))
+end
+
+discretize(chain::Chain, method::MaxLengthDiscretization) =
+  mapreduce(s -> discretize(s, method), merge, segments(chain))
+
+# -----------------
+# HELPER FUNCTIONS
+# -----------------
+
+_sides(box::Box{<:𝔼}) = sides(box)
+
+function _sides(box::Box{<:🌐})
+  A, B = extrema(box)
+  a = convert(LatLon, coords(A))
+  b = convert(LatLon, coords(B))
+  P = withcrs(box, (a.lat, b.lon), LatLon)
+
+  AP = Segment(A, P)
+  PB = Segment(P, B)
+  (_measure(AP), _measure(PB))
+end
+
+_measure(segment::Segment{<:𝔼}) = measure(segment)
+
+# TODO: Haversine returns the shortest distance between two points
+# this is not always equal to the distance between two directed points
+function _measure(segment::Segment{<:🌐})
+  T = numtype(lentype(segment))
+  🌎 = ellipsoid(datum(crs(segment)))
+  r = numconvert(T, majoraxis(🌎))
+
+  A, B = extrema(segment)
+  evaluate(Haversine(r), A, B)
+end

test/discretization.jl

@@ -425,6 +425,50 @@ end
   @test all(intersects(poly), mesh)
 end
 
+@testitem "MaxLengthDiscretization" setup = [Setup] begin
+  box = Box(cart(0, 0), cart(10, 10))
+  mesh = discretize(box, MaxLengthDiscretization(T(1)))
+  @test nvertices(mesh) == 11 * 11
+  @test nelements(mesh) == 10 * 10
+  @test eltype(mesh) <: Quadrangle
+  @test nvertices.(mesh) ⊆ [4]
+
+  box = Box(latlon(0, 0), latlon(45, 45))
+  mesh = discretize(box, MaxLengthDiscretization(T(1e5)))
+  @test nvertices(mesh) == 52 * 52
+  @test nelements(mesh) == 51 * 51
+  @test eltype(mesh) <: Quadrangle
+  @test nvertices.(mesh) ⊆ [4]
+
+  seg = Segment(cart(0, 0), cart(0, 1))
+  mesh = discretize(seg, MaxLengthDiscretization(T(0.1)))
+  @test nvertices(mesh) == 11
+  @test nelements(mesh) == 10
+  @test eltype(mesh) <: Segment
+  @test nvertices.(mesh) ⊆ [2]
+
+  seg = Segment(latlon(0, 0), latlon(0, 45))
+  mesh = discretize(seg, MaxLengthDiscretization(T(1e5)))
+  @test nvertices(mesh) == 52
+  @test nelements(mesh) == 51
+  @test eltype(mesh) <: Segment
+  @test nvertices.(mesh) ⊆ [2]
+
+  rope = Rope(cart(0, 0), cart(1, 0), cart(0, 1))
+  mesh = discretize(rope, MaxLengthDiscretization(T(0.1)))
+  @test nvertices(mesh) == 27
+  @test nelements(mesh) == 25
+  @test eltype(mesh) <: Segment
+  @test nvertices.(mesh) ⊆ [2]
+
+  ring = Ring(latlon(-45, 90), latlon(45, 90), latlon(45, -90), latlon(-45, -90))
+  mesh = discretize(ring, MaxLengthDiscretization(T(1e5)))
+  @test nvertices(mesh) == 408
+  @test nelements(mesh) == 404
+  @test eltype(mesh) <: Segment
+  @test nvertices.(mesh) ⊆ [2]
+end
+
 @testitem "Discretize" setup = [Setup] begin
   ball = Ball(cart(0, 0), T(1))
   mesh = discretize(ball)