Fix Douglas simplifies (#44)

* Update SimplifyAlgs constructors

* Start recursive DP algorithm

* Finish recursive algorithm and add test shapes

* Add prefilter argument

* Prefilter call works

* Iterative Douglas working

* Clean up and comment

* Turn on tests

* Added simplify tests

* Fix repeated last point bug

* Update comments and cleanup tests

Project.toml

@@ -21,9 +21,10 @@ ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ArchGDAL", "Distributions", "GeoFormatTypes", "GeoJSON", "LibGEOS", "Random", "Test"]
+test = ["ArchGDAL", "Distributions", "GeoFormatTypes", "GeoJSON", "JLD2", "LibGEOS", "Random", "Test"]

src/GeometryOps.jl

@@ -8,6 +8,7 @@ import Proj
 using LinearAlgebra
 import ExactPredicates
 import Proj.CoordinateTransformations.StaticArrays
+import Base.@kwdef
 
 using GeoInterface.Extents: Extents
 

src/methods/distance.jl

@@ -184,22 +184,36 @@ end
 
 # Returns the Euclidean distance between two points.
 Base.@propagate_inbounds _euclid_distance(::Type{T}, p1, p2) where T =
-    _euclid_distance(
+    sqrt(_squared_euclid_distance(T, p1, p2))
+
+# Returns the square of the euclidean distance between two points
+Base.@propagate_inbounds _squared_euclid_distance(::Type{T}, p1, p2) where T =
+    _squared_euclid_distance(
         T,
         GeoInterface.x(p1), GeoInterface.y(p1),
         GeoInterface.x(p2), GeoInterface.y(p2),
     )
 
 # Returns the Euclidean distance between two points given their x and y values.
 Base.@propagate_inbounds _euclid_distance(::Type{T}, x1, y1, x2, y2) where T =
-    T(sqrt((x2 - x1)^2 + (y2 - y1)^2))
+    sqrt(_squared_euclid_distance(T, x1, y1, x2, y2))
 
+# Returns the squared Euclidean distance between two points given their x and y values.
+Base.@propagate_inbounds _squared_euclid_distance(::Type{T}, x1, y1, x2, y2) where T =
+    T((x2 - x1)^2 + (y2 - y1)^2)
 
 #=
 Returns the minimum distance from point p0 to the line defined by endpoints p1
 and p2.
 =#
-function _distance_line(::Type{T}, p0, p1, p2) where T
+_distance_line(::Type{T}, p0, p1, p2) where T =
+    sqrt(_squared_distance_line(T, p0, p1, p2))
+
+#=
+Returns the squared minimum distance from point p0 to the line defined by
+endpoints p1 and p2.
+=#
+function _squared_distance_line(::Type{T}, p0, p1, p2) where T
     x0, y0 = GeoInterface.x(p0), GeoInterface.y(p0)
     x1, y1 = GeoInterface.x(p1), GeoInterface.y(p1)
     x2, y2 = GeoInterface.x(p2), GeoInterface.y(p2)
@@ -216,16 +230,16 @@ function _distance_line(::Type{T}, p0, p1, p2) where T
 
     c1 = sum(w .* v)
     if c1 <= 0  # p0 is closest to first endpoint
-        return _euclid_distance(T, x0, y0, xfirst, yfirst)
+        return _squared_euclid_distance(T, x0, y0, xfirst, yfirst)
     end
 
     c2 = sum(v .* v)
     if c2 <= c1 # p0 is closest to last endpoint
-        return _euclid_distance(T, x0, y0, xlast, ylast)
+        return _squared_euclid_distance(T, x0, y0, xlast, ylast)
     end
 
     b2 = c1 / c2  # projection fraction
-    return _euclid_distance(T, x0, y0, xfirst + (b2 * v[1]), yfirst + (b2 * v[2]))
+    return _squared_euclid_distance(T, x0, y0, xfirst + (b2 * v[1]), yfirst + (b2 * v[2]))
 end