diff --git a/src/Meshes.jl b/src/Meshes.jl index 02514fe24..fb3126c55 100644 --- a/src/Meshes.jl +++ b/src/Meshes.jl @@ -15,9 +15,8 @@ using Random using Bessels: gamma using Unitful: AbstractQuantity, numtype using StatsBase: AbstractWeights, Weights, quantile -using Distances: PreMetric, Euclidean, Mahalanobis +using Distances: PreMetric, MinkowskiMetric, Euclidean, Mahalanobis using Distances: Haversine, SphericalAngle -using Distances: evaluate, result_type using Rotations: Rotation, QuatRotation, Angle2d using Rotations: rotation_between using CoordRefSystems: Basic, Projected, Geographic @@ -37,8 +36,8 @@ import Base: ==, ! import Base: +, -, * import Base: <, >, ≤, ≥ import StatsBase: sample -import Distances: evaluate -import NearestNeighbors: MinkowskiMetric +import Distances: evaluate, result_type +using NearestNeighbors: eval_pow, eval_diff # Transforms API import TransformsBase: Transform, → @@ -558,6 +557,7 @@ export measurematrix, adjacencymatrix, atol, + EuclideanDistance, # visualization viz, diff --git a/src/distances.jl b/src/distances.jl index 569f063c6..eb1618c34 100644 --- a/src/distances.jl +++ b/src/distances.jl @@ -5,6 +5,51 @@ # flip arguments so that points always come first evaluate(d::PreMetric, g::Geometry, p::Point) = evaluate(d, p, g) +# ---------- +# EUCLIDEAN +# ---------- + +struct EuclideanDistance <: MinkowskiMetric end + +# Distances.jl interface +result_type(::EuclideanDistance, ::Type{T₁}, ::Type{T₂}) where {T₁,T₂} = float(promote_type(T₁, T₂)) + +@propagate_inbounds function (d::EuclideanDistance)(a, b) + @boundscheck if length(a) ≠ length(b) + throw( + DimensionMismatch( + "first collection has length $(length(a)) which does not match the length of the second, $(length(b))." + ) + ) + end + norm((bᵢ - aᵢ) for (aᵢ, bᵢ) in zip(a, b)) +end + +# implementations for geometries +(d::EuclideanDistance)(p₁::Point, p₂::Point) = norm(p₂ - p₁) + +function (d::EuclideanDistance)(p::Point, l::Line) + a, b = l(0), l(1) + u = p - a + v = b - a + α = (u ⋅ v) / (v ⋅ v) + norm(u - α * v) +end + +(d::EuclideanDistance)(l::Line, p::Point) = d(p, l) + +function (d::EuclideanDistance)(l₁::Line, l₂::Line) + λ₁, λ₂, r, rₐ = intersectparameters(l₁(0), l₁(1), l₂(0), l₂(1)) + + if (r == rₐ == 2) || (r == rₐ == 1) # lines intersect or are colinear + zero(result_type(d, lentype(l₁), lentype(l₂))) + elseif (r == 1) && (rₐ == 2) # lines are parallel + d(l₁(0), l₂) + else # get distance between closest points on each line + d(l₁(λ₁), l₂(λ₂)) + end +end + """ evaluate(distance::Euclidean, point, line)