|
| 1 | +#= |
| 2 | +# Smooth |
| 3 | +
|
| 4 | +Geometry smoothing is meant to make shapes more aesthetically pleasing, usually by rounding out rough edges and corners. |
| 5 | +
|
| 6 | +You can do this by the [`smooth`](@ref) function, which uses the [`Chaikin`](@ref) algorithm by default. |
| 7 | +
|
| 8 | +## Example |
| 9 | +
|
| 10 | +```@example smooth |
| 11 | +using CairoMakie |
| 12 | +import GeoInterface as GI, GeometryOps as GO |
| 13 | +
|
| 14 | +line = GI.LineString([(0.0, 0.0), (1.0, 1.0), (2.0, 0.0)]) |
| 15 | +smoothed = GO.smooth(line) |
| 16 | +smoothed_2 = GO.smooth(line; iterations=2) |
| 17 | +
|
| 18 | +f, a, p = lines(line; label = "Original") |
| 19 | +lines!(a, smoothed; label = "1 iteration") |
| 20 | +lines!(a, smoothed_2; label = "2 iterations") |
| 21 | +axislegend(a) |
| 22 | +fig |
| 23 | +``` |
| 24 | +
|
| 25 | +Smoothing also works on the [`Spherical`](@ref) manifold, |
| 26 | +similarly to the planar manifold (default): |
| 27 | +```@example smooth |
| 28 | +using CairoMakie |
| 29 | +import GeoInterface as GI, GeometryOps as GO |
| 30 | +
|
| 31 | +line = GI.LineString([(0.0, 0.0), (1.0, 1.0), (2.0, 0.0)]) |
| 32 | +smoothed = GO.smooth(GO.Spherical(), line) |> x -> GO.transform(GO.UnitSpherical.GeographicFromUnitSpherical(), x) |
| 33 | +smoothed_2 = GO.smooth(GO.Spherical(), line; iterations=2) |> x -> GO.transform(GO.UnitSpherical.GeographicFromUnitSpherical(), x) |
| 34 | +
|
| 35 | +f, a, p = lines(line; label = "Original", axis = (; title = "Spherical smoothing")) |
| 36 | +lines!(a, smoothed; label = "1 iteration") |
| 37 | +lines!(a, smoothed_2; label = "2 iterations") |
| 38 | +axislegend(a) |
| 39 | +fig |
| 40 | +``` |
| 41 | +=# |
| 42 | +""" |
| 43 | + Chaikin(; iterations=1, manifold=Planar()) |
| 44 | +
|
| 45 | +Smooths geometries using Chaikin's corner-cutting algorithm [^1]. |
| 46 | +This algorithm "slices" off every corner of the geometry to smooth it out, |
| 47 | +equivalent to a sequence of quadratic Bezier curves. |
| 48 | +
|
| 49 | +## Keywords |
| 50 | +- `iterations`: the number of times to apply the algorithm. |
| 51 | +- `manifold`: the [`Manifold`](@ref) to smooth the geometry on. Currently, [`Planar`](@ref) and [`Spherical`](@ref) are supported. |
| 52 | +
|
| 53 | +# Extended help |
| 54 | +
|
| 55 | +The algorithm is very simple; for each corner of the line (a -> b -> c), |
| 56 | +insert two new points and remove b, such that `a -> b -> c` becomes |
| 57 | +`a -> q -> r -> c`, where `q` and `r` are the new points such that: |
| 58 | +
|
| 59 | +```math |
| 60 | +q = 3/4 * b + 1/4 * a |
| 61 | +r = 3/4 * b + 1/4 * c |
| 62 | +``` |
| 63 | +
|
| 64 | +In practice the replacement happens on the level of each edge. |
| 65 | +
|
| 66 | +## References |
| 67 | +[^1]: Chaikin, G. An algorithm for high speed curve generation. Computer Graphics and Image Processing 3 (1974), 346-349 |
| 68 | +""" |
| 69 | +@kwdef struct Chaikin{M} <: Algorithm{M} |
| 70 | + manifold::M = Planar() |
| 71 | + iterations::Int = 1 |
| 72 | +end |
| 73 | + |
| 74 | +""" |
| 75 | + smooth(alg::Algorithm, geom) |
| 76 | + smooth(geom; kw...) |
| 77 | +
|
| 78 | +Smooths a geometry using the provided algorithm. |
| 79 | +
|
| 80 | +The default algorithm is [`Chaikin()`](@ref), which can be used on the spherical or planar manifolds. |
| 81 | +""" |
| 82 | +smooth(geom; kw...) = smooth(Chaikin(; kw...), geom) |
| 83 | +smooth(m::Manifold, geom; kw...) = smooth(Chaikin(; manifold=m, kw...), geom) |
| 84 | +function smooth(alg::Algorithm, geom; kw...) |
| 85 | + _smooth_function(trait, geom) = _smooth(alg, trait, geom) |
| 86 | + return apply( |
| 87 | + WithTrait(_smooth_function), |
| 88 | + TraitTarget{Union{GI.AbstractCurveTrait,GI.MultiPointTrait,GI.PointTrait}}(), |
| 89 | + geom; |
| 90 | + kw... |
| 91 | + ) |
| 92 | +end |
| 93 | + |
| 94 | +_smooth(alg, ::GI.PointTrait, geom) = geom |
| 95 | +_smooth(alg, ::GI.MultiPointTrait, geom) = geom |
| 96 | + |
| 97 | +function _smooth(alg::Chaikin{<: Planar}, trait::Trait, geom) where {M, Trait <: Union{GI.LineStringTrait,GI.LinearRingTrait}} |
| 98 | + isring = Trait <: GI.LinearRingTrait |
| 99 | + points = tuple_points(geom) |
| 100 | + if isring && first(points) != last(points) |
| 101 | + push!(points, first(points)) |
| 102 | + end |
| 103 | + smoothed_points = _chaikin_smooth(alg.manifold, points, alg.iterations, isring) |
| 104 | + return rebuild(geom, smoothed_points) |
| 105 | +end |
| 106 | + |
| 107 | +function _smooth(alg::Chaikin{<: M}, trait::Trait, geom) where {M <: Spherical, Trait <: Union{GI.LineStringTrait,GI.LinearRingTrait}} |
| 108 | + isring = Trait <: GI.LinearRingTrait |
| 109 | + points = apply(UnitSphereFromGeographic(), GI.PointTrait(), geom).geom |
| 110 | + if isring && first(points) != last(points) |
| 111 | + push!(points, first(points)) |
| 112 | + end |
| 113 | + smoothed_points = _chaikin_smooth(alg.manifold, points, alg.iterations, isring) |
| 114 | + return rebuild(geom, smoothed_points) |
| 115 | +end |
| 116 | + |
| 117 | +function _chaikin_smooth(manifold::M, points::Vector{P}, iterations::Int, isring::Bool) where {M <: Manifold, P} |
| 118 | + # points is expected to be a vector of points |
| 119 | + smoothed_points = points |
| 120 | + for itr in 1:iterations |
| 121 | + num_points = length(smoothed_points) |
| 122 | + if isring |
| 123 | + n = 1 |
| 124 | + new_points = Vector{P}(undef, num_points * 2 - 1) |
| 125 | + else |
| 126 | + n = 2 |
| 127 | + # Need to add the first point |
| 128 | + new_points = Vector{P}(undef, num_points * 2) |
| 129 | + new_points[begin] = smoothed_points[begin] |
| 130 | + new_points[end] = smoothed_points[end] |
| 131 | + end |
| 132 | + # fill!(new_points, (P <: NTuple{2, Float64} ? (-9999.0, -9999.0) : UnitSphericalPoint(-9999.0, -9999.0, -9999.0))) |
| 133 | + |
| 134 | + for i in eachindex(smoothed_points)[begin:end-1] |
| 135 | + p1 = smoothed_points[i] |
| 136 | + p2 = smoothed_points[i+1] |
| 137 | + _add_smoothed_points!(manifold, new_points, p1, p2, n) |
| 138 | + n += 2 |
| 139 | + end |
| 140 | + |
| 141 | + if isring # Close it |
| 142 | + new_points[end] = new_points[begin] |
| 143 | + end |
| 144 | + |
| 145 | + smoothed_points = new_points |
| 146 | + end |
| 147 | + |
| 148 | + return smoothed_points |
| 149 | +end |
| 150 | + |
| 151 | +function _add_smoothed_points!(::Planar, new_points, p1, p2, n) |
| 152 | + q_x = 0.75 * GI.x(p1) + 0.25 * GI.x(p2) |
| 153 | + q_y = 0.75 * GI.y(p1) + 0.25 * GI.y(p2) |
| 154 | + r_x = 0.25 * GI.x(p1) + 0.75 * GI.x(p2) |
| 155 | + r_y = 0.25 * GI.y(p1) + 0.75 * GI.y(p2) |
| 156 | + |
| 157 | + new_points[n] = (q_x, q_y) |
| 158 | + new_points[n+1] = (r_x, r_y) |
| 159 | +end |
| 160 | + |
| 161 | +# For spherical points, we can simply slerp. |
| 162 | +function _add_smoothed_points!(::Spherical, new_points, p1, p2, n) |
| 163 | + q = slerp(p1, p2, 0.25) |
| 164 | + r = slerp(p1, p2, 0.75) |
| 165 | + |
| 166 | + new_points[n] = q |
| 167 | + new_points[n+1] = r |
| 168 | +end |
0 commit comments