Refactor simplification methods (#1019)

* Initial refactor

* Add BinarySearchSimplification

* More adjustments

* Fix docstring

* More adjustments

* More adjustments

* More adjustments

* More adjustments

* More adjustments

* Formatting

* Minor adjustments

Author: juliohm

docs/src/algorithms/simplification.md

@@ -7,86 +7,103 @@ import CairoMakie as Mke # hide
 
 ```@docs
 simplify
-decimate
 SimplificationMethod
 ```
 
-## Douglas-Peucker
+## SelingerSimplification
 
 ```@docs
-DouglasPeucker
+SelingerSimplification
 ```
 
 ```@example simplification
-# polygonal area
-polyarea = PolyArea([(0.22926679, 0.47329807), (0.23094065, 0.44913536), (0.2569517, 0.38217533),
-                     (0.3072999, 0.272418), (0.34814754, 0.18421611), (0.37949452, 0.11756973),
-                     (0.4013409, 0.07247882), (0.41368666, 0.048943404), (0.42597583, 0.031655528),
-                     (0.4382084, 0.0206152), (0.45038435, 0.015822414), (0.4625037, 0.017277176),
-                     (0.47175184, 0.02439156), (0.47812873, 0.03716557), (0.4816344, 0.055599205),
-                     (0.48226887, 0.07969247), (0.48172843, 0.10446181), (0.4800131, 0.12990724),
-                     (0.47712287, 0.15602873), (0.47305775, 0.18282633), (0.47093934, 0.20558843),
-                     (0.47076762, 0.22431506), (0.47254258, 0.23900622), (0.47626427, 0.24966191),
-                     (0.47768936, 0.25845313), (0.47681788, 0.26537988), (0.4736498, 0.27044216),
-                     (0.46818516, 0.27363995), (0.4613889, 0.27232954), (0.45326096, 0.2665109),
-                     (0.44380143, 0.256184), (0.43301025, 0.24134888), (0.4246466, 0.22978415),
-                     (0.41871038, 0.22148979), (0.4152017, 0.21646582), (0.4141205, 0.21471222),
-                     (0.41227448, 0.21589448), (0.40966362, 0.22001258), (0.40628797, 0.22706655),
-                     (0.40214747, 0.23705636), (0.40200475, 0.24653101), (0.40585983, 0.25549048),
-                     (0.41371268, 0.2639348), (0.4255633, 0.2718639), (0.4378565, 0.28495985),
-                     (0.4505922, 0.30322257), (0.46377045, 0.32665208), (0.47739124, 0.35524836),
-                     (0.5046394, 0.36442512), (0.5455148, 0.35418236), (0.60001767, 0.32452005),
-                     (0.66814786, 0.27543822), (0.7186763, 0.24664374), (0.75160307, 0.23813659),
-                     (0.76692814, 0.2499168), (0.7646515, 0.28198436), (0.7769703, 0.29925033),
-                     (0.8038847, 0.3017147), (0.84539455, 0.28937748), (0.9015, 0.26223865),
-                     (0.94408435, 0.24899776), (0.9731477, 0.24965483), (0.98869, 0.26420987),
-                     (0.9907113, 0.29266283), (0.9849871, 0.31338844), (0.97151726, 0.32638666),
-                     (0.950302, 0.3316575), (0.9213412, 0.32920095), (0.8798396, 0.34078467),
-                     (0.8257972, 0.36640862), (0.7592141, 0.40607283), (0.6800901, 0.4597773),
-                     (0.6450007, 0.49104902), (0.6539457, 0.49988794), (0.7069251, 0.48629412),
-                     (0.803939, 0.45026752), (0.877913, 0.4226481), (0.9288472, 0.40343583),
-                     (0.9567415, 0.39263073), (0.961596, 0.39023277), (0.9419039, 0.40523484),
-                     (0.89766514, 0.43763688), (0.8288798, 0.48743892), (0.7355478, 0.55464095),
-                     (0.6655121, 0.60063523), (0.6187727, 0.6254217), (0.5953296, 0.62900037),
-                     (0.5951828, 0.6113712), (0.57516366, 0.60261106), (0.53527224, 0.6027198),
-                     (0.4755085, 0.6116975), (0.3958725, 0.6295441), (0.33913234, 0.6398651),
-                     (0.30528808, 0.6426605), (0.2943397, 0.6379303), (0.30628717, 0.6256744),
-                     (0.32149008, 0.6093727), (0.33994842, 0.5890249), (0.36166218, 0.5646312),
-                     (0.38663134, 0.5361916), (0.3919681, 0.520893), (0.3776725, 0.5187355),
-                     (0.34374446, 0.52971905), (0.29018405, 0.5538437), (0.25439468, 0.5678829),
-                     (0.2363764, 0.5718367), (0.23612918, 0.56570506), (0.25365302, 0.549488),
-                     (0.2733971, 0.5246488), (0.29536137, 0.49118724), (0.3195459, 0.4491035),
-                     (0.34595063, 0.39839754), (0.3647463, 0.3590396), (0.37593287, 0.33102974),
-                     (0.37951034, 0.31436795), (0.37547874, 0.30905423), (0.36070493, 0.3204269),
-                     (0.33518887, 0.348486), (0.29893062, 0.3932315), (0.25193012, 0.45466346)])
+poly = PolyArea([(0.22926679, 0.47329807), (0.23094065, 0.44913536), (0.2569517, 0.38217533),
+                 (0.3072999, 0.272418), (0.34814754, 0.18421611), (0.37949452, 0.11756973),
+                 (0.4013409, 0.07247882), (0.41368666, 0.048943404), (0.42597583, 0.031655528),
+                 (0.4382084, 0.0206152), (0.45038435, 0.015822414), (0.4625037, 0.017277176),
+                 (0.47175184, 0.02439156), (0.47812873, 0.03716557), (0.4816344, 0.055599205),
+                 (0.48226887, 0.07969247), (0.48172843, 0.10446181), (0.4800131, 0.12990724),
+                 (0.47712287, 0.15602873), (0.47305775, 0.18282633), (0.47093934, 0.20558843),
+                 (0.47076762, 0.22431506), (0.47254258, 0.23900622), (0.47626427, 0.24966191),
+                 (0.47768936, 0.25845313), (0.47681788, 0.26537988), (0.4736498, 0.27044216),
+                 (0.46818516, 0.27363995), (0.4613889, 0.27232954), (0.45326096, 0.2665109),
+                 (0.44380143, 0.256184), (0.43301025, 0.24134888), (0.4246466, 0.22978415),
+                 (0.41871038, 0.22148979), (0.4152017, 0.21646582), (0.4141205, 0.21471222),
+                 (0.41227448, 0.21589448), (0.40966362, 0.22001258), (0.40628797, 0.22706655),
+                 (0.40214747, 0.23705636), (0.40200475, 0.24653101), (0.40585983, 0.25549048),
+                 (0.41371268, 0.2639348), (0.4255633, 0.2718639), (0.4378565, 0.28495985),
+                 (0.4505922, 0.30322257), (0.46377045, 0.32665208), (0.47739124, 0.35524836),
+                 (0.5046394, 0.36442512), (0.5455148, 0.35418236), (0.60001767, 0.32452005),
+                 (0.66814786, 0.27543822), (0.7186763, 0.24664374), (0.75160307, 0.23813659),
+                 (0.76692814, 0.2499168), (0.7646515, 0.28198436), (0.7769703, 0.29925033),
+                 (0.8038847, 0.3017147), (0.84539455, 0.28937748), (0.9015, 0.26223865),
+                 (0.94408435, 0.24899776), (0.9731477, 0.24965483), (0.98869, 0.26420987),
+                 (0.9907113, 0.29266283), (0.9849871, 0.31338844), (0.97151726, 0.32638666),
+                 (0.950302, 0.3316575), (0.9213412, 0.32920095), (0.8798396, 0.34078467),
+                 (0.8257972, 0.36640862), (0.7592141, 0.40607283), (0.6800901, 0.4597773),
+                 (0.6450007, 0.49104902), (0.6539457, 0.49988794), (0.7069251, 0.48629412),
+                 (0.803939, 0.45026752), (0.877913, 0.4226481), (0.9288472, 0.40343583),
+                 (0.9567415, 0.39263073), (0.961596, 0.39023277), (0.9419039, 0.40523484),
+                 (0.89766514, 0.43763688), (0.8288798, 0.48743892), (0.7355478, 0.55464095),
+                 (0.6655121, 0.60063523), (0.6187727, 0.6254217), (0.5953296, 0.62900037),
+                 (0.5951828, 0.6113712), (0.57516366, 0.60261106), (0.53527224, 0.6027198),
+                 (0.4755085, 0.6116975), (0.3958725, 0.6295441), (0.33913234, 0.6398651),
+                 (0.30528808, 0.6426605), (0.2943397, 0.6379303), (0.30628717, 0.6256744),
+                 (0.32149008, 0.6093727), (0.33994842, 0.5890249), (0.36166218, 0.5646312),
+                 (0.38663134, 0.5361916), (0.3919681, 0.520893), (0.3776725, 0.5187355),
+                 (0.34374446, 0.52971905), (0.29018405, 0.5538437), (0.25439468, 0.5678829),
+                 (0.2363764, 0.5718367), (0.23612918, 0.56570506), (0.25365302, 0.549488),
+                 (0.2733971, 0.5246488), (0.29536137, 0.49118724), (0.3195459, 0.4491035),
+                 (0.34595063, 0.39839754), (0.3647463, 0.3590396), (0.37593287, 0.33102974),
+                 (0.37951034, 0.31436795), (0.37547874, 0.30905423), (0.36070493, 0.3204269),
+                 (0.33518887, 0.348486), (0.29893062, 0.3932315), (0.25193012, 0.45466346)])
 
-simp1 = simplify(polyarea, DouglasPeucker(0.01))
-simp2 = simplify(polyarea, DouglasPeucker(0.05))
-simp3 = simplify(polyarea, DouglasPeucker(0.10))
+simp1 = simplify(poly, SelingerSimplification(0.01))
+simp2 = simplify(poly, SelingerSimplification(0.05))
+simp3 = simplify(poly, SelingerSimplification(0.10))
 
 fig = Mke.Figure(size = (800, 800))
-viz(fig[1,1], polyarea)
+viz(fig[1,1], poly)
 viz(fig[1,2], simp1)
 viz(fig[2,1], simp2)
 viz(fig[2,2], simp3)
 fig
 ```
 
-## Selinger
+## DouglasPeuckerSimplification
 
 ```@docs
-Selinger
+DouglasPeuckerSimplification
 ```
 
 ```@example simplification
-simp1 = simplify(polyarea, Selinger(0.01))
-simp2 = simplify(polyarea, Selinger(0.05))
-simp3 = simplify(polyarea, Selinger(0.10))
+simp1 = simplify(poly, DouglasPeuckerSimplification(0.01))
+simp2 = simplify(poly, DouglasPeuckerSimplification(0.05))
+simp3 = simplify(poly, DouglasPeuckerSimplification(0.10))
 
 fig = Mke.Figure(size = (800, 800))
-viz(fig[1,1], polyarea)
+viz(fig[1,1], poly)
 viz(fig[1,2], simp1)
 viz(fig[2,1], simp2)
 viz(fig[2,2], simp3)
 fig
-```
+```
+
+## MinMaxSimplification
+
+```@docs
+MinMaxSimplification
+```
+
+```@example simplification
+simp1 = simplify(poly, MinMaxSimplification(DouglasPeuckerSimplification, max=20))
+simp2 = simplify(poly, MinMaxSimplification(DouglasPeuckerSimplification, max=10))
+simp3 = simplify(poly, MinMaxSimplification(DouglasPeuckerSimplification, max=5))
+
+fig = Mke.Figure(size = (800, 800))
+viz(fig[1,1], poly)
+viz(fig[1,2], simp1)
+viz(fig[2,1], simp2)
+viz(fig[2,2], simp3)
+fig
+```

src/Meshes.jl

@@ -457,10 +457,10 @@ export
 
   # simplification
   SimplificationMethod,
-  DouglasPeucker,
-  Selinger,
+  SelingerSimplification,
+  DouglasPeuckerSimplification,
+  MinMaxSimplification,
   simplify,
-  decimate,
 
   # bounding boxes
   boundingbox,

src/simplification.jl

@@ -12,15 +12,11 @@ abstract type SimplificationMethod end
 """
     simplify(object, method)
 
-Simplify `object` with given `method`.
-
-See also [`decimate`](@ref).
+Simplify geometric `object` with given `method`.
 """
 function simplify end
 
-simplify(box::Box, method::SimplificationMethod) = _simplify(box, Val(embeddim(box)), method)
-
-_simplify(box::Box, ::Val{2}, method::SimplificationMethod) = PolyArea(simplify(boundary(box), method))
+simplify(box::Box{𝔼{2}}, method::SimplificationMethod) = PolyArea(simplify(boundary(box), method))
 
 simplify(polygon::Polygon, method::SimplificationMethod) = PolyArea([simplify(ring, method) for ring in rings(polygon)])
 
@@ -32,22 +28,6 @@ simplify(domain::Domain, method::SimplificationMethod) = GeometrySet([simplify(e
 # IMPLEMENTATIONS
 # ----------------
 
-include("simplification/douglaspeucker.jl")
 include("simplification/selinger.jl")
-
-# ----------
-# UTILITIES
-# ----------
-
-"""
-    decimate(object, [ϵ]; min=3, max=typemax(Int), maxiter=10)
-
-Simplify `object` with an appropriate simplification method
-and deviation tolerance `ϵ`.
-
-If the tolerance `ϵ` is not provided, perform binary search until
-the number of vertices is between `min` and `max` or until the
-number of iterations reaches a maximum `maxiter`.
-"""
-decimate(object, ϵ=nothing; min=3, max=typemax(Int), maxiter=10) =
-  simplify(object, DouglasPeucker(ϵ, min=min, max=max, maxiter=maxiter))
+include("simplification/douglaspeucker.jl")
+include("simplification/minmax.jl")

src/simplification/douglaspeucker.jl

@@ -3,84 +3,34 @@
 # ------------------------------------------------------------------
 
 """
-    DouglasPeucker([ϵ]; min=3, max=typemax(Int), maxiter=10)
+    DouglasPeuckerSimplification(τ)
 
-Simplify geometries with Douglas-Peucker algorithm. The higher
-is the tolerance `ϵ`, the more aggressive is the simplification.
+Douglas-Peucker's simplification algorithm with tolerance `τ` in length units
+(default to meter).
 
-If the tolerance `ϵ` is not provided, perform binary search until
-the number of vertices is between `min` and `max` or until the
-number of iterations reaches a maximum `maxiter`.
+The higher is the tolerance, the more aggressive is the simplification.
 
 ## References
 
 * Douglas, D. and Peucker, T. 1973. [Algorithms for the Reduction of
   the Number of Points Required to Represent a Digitized Line or its
   Caricature](https://www.sciencedirect.com/science/article/abs/pii/0167839691900198)
 """
-struct DouglasPeucker{T} <: SimplificationMethod
-  ϵ::T
-  min::Int
-  max::Int
-  maxiter::Int
+struct DouglasPeuckerSimplification{ℒ<:Len} <: SimplificationMethod
+  τ::ℒ
+  DouglasPeuckerSimplification(τ::ℒ) where {ℒ<:Len} = new{float(ℒ)}(τ)
 end
 
-DouglasPeucker(ϵ=nothing; min=3, max=typemax(Int), maxiter=10) = DouglasPeucker(_ϵ(ϵ), min, max, maxiter)
+DouglasPeuckerSimplification(τ) = DouglasPeuckerSimplification(addunit(τ, u"m"))
 
-_ϵ(ϵ::Nothing) = ϵ
-_ϵ(ϵ::Number) = _ϵ(addunit(ϵ, u"m"))
-_ϵ(ϵ::Len) = float(ϵ)
-
-function simplify(chain::Chain, method::DouglasPeucker)
-  v = if isnothing(method.ϵ)
-    # perform binary search with other parameters
-    βsimplify(vertices(chain), method.min, method.max, method.maxiter)
-  else
-    # perform Douglas-Peucker ϵ-simplification
-    ϵsimplify(vertices(chain), method.ϵ)
-  end |> collect
-  isclosed(chain) ? Ring(v) : Rope(v)
-end
-
-# simplification by means of binary search
-function βsimplify(v::AbstractVector{P}, min, max, maxiter) where {P<:Point}
-  i = 0
-  u = v
-  n = length(u)
-  a = zero(lentype(P))
-  b = initeps(u)
-  while !(min ≤ n ≤ max) && i < maxiter
-    # midpoint candidate
-    ϵ = (a + b) / 2
-
-    # evaluate at midpoint
-    u = ϵsimplify(v, ϵ)
-    n = length(u)
-
-    # binary search
-    n < min && (b = ϵ)
-    n > max && (a = ϵ)
-
-    i += 1
-  end
-
-  u
-end
-
-# initial ϵ guess for a given chain
-function initeps(v::AbstractVector{P}) where {P<:Point}
-  n = length(v)
-  ϵ = typemax(lentype(P))
-  l = Line(first(v), last(v))
-  d = [evaluate(Euclidean(), v[i], l) for i in 2:(n - 1)]
-  ϵ = quantile(d, 0.25)
-  2ϵ
+function simplify(chain::Chain, method::DouglasPeuckerSimplification)
+  verts = _douglaspeucker(vertices(chain), method.τ) |> collect
+  isclosed(chain) ? Ring(verts) : Rope(verts)
 end
 
 # simplify chain assuming it is open
-function ϵsimplify(v::AbstractVector{P}, ϵ) where {P<:Point}
-  # find vertex with maximum distance
-  # to reference line
+function _douglaspeucker(v::AbstractVector{P}, τ) where {P<:Point}
+  # find vertex with maximum distance to reference line
   l = Line(first(v), last(v))
   imax, dmax = 0, zero(lentype(P))
   for i in 2:(length(v) - 1)
@@ -91,11 +41,11 @@ function ϵsimplify(v::AbstractVector{P}, ϵ) where {P<:Point}
     end
   end
 
-  if dmax < ϵ
+  if dmax < τ
     [first(v), last(v)]
   else
-    v₁ = ϵsimplify(v[begin:imax], ϵ)
-    v₂ = ϵsimplify(v[imax:end], ϵ)
+    v₁ = _douglaspeucker(v[begin:imax], τ)
+    v₂ = _douglaspeucker(v[imax:end], τ)
     [v₁[begin:(end - 1)]; v₂]
   end
 end

src/simplification/minmax.jl

@@ -0,0 +1,54 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    MinMaxSimplification(method; min=3, max=typemax(Int), maxiter=10)
+
+Simplify geometries with binary search algorithm and a parent simplification `method`.
+
+The simplification is performed until the number of vertices is in the `[min, max]`
+range or until a maximum number of iterations `maxiter` is reached.
+"""
+struct MinMaxSimplification{M} <: SimplificationMethod
+  method::M
+  min::Int
+  max::Int
+  maxiter::Int
+end
+
+MinMaxSimplification(method; min=3, max=typemax(Int), maxiter=10) = MinMaxSimplification(method, min, max, maxiter)
+
+function simplify(c::Chain, m::MinMaxSimplification)
+  i = 0
+  s = c
+  n = nvertices(c)
+  a, b = _initrange(c)
+  while !(m.min ≤ n ≤ m.max) && i < m.maxiter
+    # midpoint candidate
+    τ = (a + b) / 2
+
+    # evaluate at midpoint
+    s = simplify(c, m.method(τ))
+    n = nvertices(s)
+
+    # binary search
+    n < m.min && (b = τ)
+    n > m.max && (a = τ)
+
+    i += 1
+  end
+
+  s
+end
+
+# initial range for binary search
+function _initrange(c)
+  v = vertices(c)
+  n = length(v)
+  l = Line(first(v), last(v))
+  d = [evaluate(Euclidean(), v[i], l) for i in 2:(n - 1)]
+  z = zero(lentype(c))
+  τ = quantile(d, 0.25)
+  (z, 2τ)
+end