rework continuous frechet distance lower bound and fix simplification nomenclature

Author: derohde

README.md

@@ -25,8 +25,8 @@ By default, Fred will automatically determine the number of threads to use. If y
 - returns: `fred.Continuous_Frechet_Result` with members `value`, `time_bounds`: running-time for upper and lower bound, `number_searches`: number of free space diagrams built, `time_searches`: running-time for free spaces
 
 ###### continuous Frechet distance config
-- approximation error: `fred.set_continuous_frechet_epsilon(epsilon)` with parameter `epsilon`, which defaults to 0.001
-- rounding (rounding up to 3 decimals): `fred.set_continuous_frechet_rounding(round)` with parameter `round`, which defaults to true
+- approximation error in percent of distance: `fred.set_continuous_frechet_error(double percent)` with parameter `percent`, which defaults to 1
+- rounding: `fred.set_continuous_frechet_rounding(round)` with parameter `round`, which defaults to true
 
 #### discrete Fréchet distance
 - signature: `fred.discrete_frechet(curve1, curve2)`
@@ -40,19 +40,19 @@ By default, Fred will automatically determine the number of threads to use. If y
 
 All simplifications are vertex-restricted!
 
-#### weak minimum error simplification
+#### minimum error simplification
 - graph approach from [**Polygonal Approximations of a Curve — Formulations and Algorithms**](https://www.sciencedirect.com/science/article/pii/B9780444704672500114)
-- signature: `fred.weak_minimum_error_simplification(fred.Curve, int complexity)`
+- signature: `fred.minimum_error_simplification(fred.Curve, int complexity)`
 - returns: `fred.Curve`that uses input curves vertices, with `complexity` number of vertices and that has minimum distance to input curve
 
-#### approximate weak minimum link simplification
+#### approximate minimum link simplification
 - algorithm "FS" from [**Near-Linear Time Approximation Algorithms for Curve Simplification**](https://link.springer.com/article/10.1007/s00453-005-1165-y)
-- signature: `fred.approximate_weak_minimum_error_simplification(fred.Curve, double error)`
+- signature: `fred.approximate_minimum_link_simplification(fred.Curve, double error)`
 - returns: `fred.Curve` that uses input curves vertices, is of small complexity and with distance to input curve at most `error`
 
-#### approximate weak minimum error simplification
-- binary search on `fred.approximate_weak_minimum_link_simplification`
-- signature: `fred.approximate_weak_minimum_error_simplification(fred.Curve, int complexity)`
+#### approximate minimum error simplification
+- binary search on `fred.approximate_minimum_link_simplification`
+- signature: `fred.approximate_minimum_error_simplification(fred.Curve, int complexity)`
 - returns: `fred.Curve`that uses input curves vertices, with `complexity` number of vertices and that has small distance to input curve
 
 ### Clustering
@@ -68,7 +68,7 @@ A `fred.Distance_Matrix()` can be used to speed up consecutive calls of `fred.di
     - `l`: maximum complexity of the centers
     - `distances`: `fred.Distance_Matrix`, defaults to empty `fred.Distance_Matrix`
     - `random_first_center`: determines if first center is chosen uniformly at random or first curve is used as first center, optional, defaults to true
-    - `fast_simplification`: determines whether to use the weak minimum error simplification or the faster approximate weak minimum error simplification, defaults to false
+    - `fast_simplification`: determines whether to use the minimum error simplification or the faster approximate minimum error simplification, defaults to false
 - returns: `fred.Clustering_Result` with mebers 
     - `value`: objective value 
     - `time`: running-time 
@@ -80,7 +80,7 @@ A `fred.Distance_Matrix()` can be used to speed up consecutive calls of `fred.di
     - `k`: number of centers
     - `l`: maximum complexity of the centers
     - `distances`: `fred.Distance_Matrix`, defaults to empty `fred.Distance_Matrix`
-    - `fast_simplification`: determines whether to use the weak minimum error simplification or the faster approximate weak minimum error simplification, defaults to false
+    - `fast_simplification`: determines whether to use the minimum error simplification or the faster approximate minimum error simplification, defaults to false
 - returns: `fred.Clustering_Result` with mebers 
     - `value`: objective value 
     - `time`: running-time 
@@ -136,7 +136,7 @@ curve2d2 = fred.Curve(np.array([[1., -1.], [2., -2.], [3., -1.]]), "optional nam
 print(curve2d1)
 
 Fred.plot_curve(curve2d1, curve2d2)
-Fred.plot_curve(curve2d2, fred.weak_minimum_error_simplification(curve2d2, 2))
+Fred.plot_curve(curve2d2, fred.minimum_error_simplification(curve2d2, 2))
 
 print("distance is {}".format(fred.continuous_frechet(curve2d1, curve2d2).value))
 

include/clustering.hpp

@@ -140,12 +140,12 @@ Clustering_Result kl_cluster(const curve_number_t num_centers, const curve_size_
 
     auto simplify = [&](const curve_number_t i) {
         if (fast_simplification) {
-            auto simplified_curve = Simplification::approximate_weak_minimum_error_simplification(const_cast<Curve&>(in[i]), ell);
+            auto simplified_curve = Simplification::approximate_minimum_error_simplification(const_cast<Curve&>(in[i]), ell);
             simplified_curve.set_name("Simplification of " + in[i].get_name());
             return simplified_curve;
         } else {
             Simplification::Subcurve_Shortcut_Graph graph(const_cast<Curve&>(in[i]));
-            auto simplified_curve = graph.weak_minimum_error_simplification(ell);
+            auto simplified_curve = graph.minimum_error_simplification(ell);
             simplified_curve.set_name("Simplification of " + in[i].get_name());
             return simplified_curve;
         }
@@ -268,7 +268,7 @@ Clustering_Result one_median_sampling(const curve_size_t ell, const Curves &in,
 
         for (curve_number_t i = 0; i < in.size(); ++i) {
             Simplification::Subcurve_Shortcut_Graph graph(const_cast<Curve&>(in[i]));
-            auto simplified_curve = graph.weak_minimum_error_simplification(ell);
+            auto simplified_curve = graph.minimum_error_simplification(ell);
             simplified_curve.set_name("Simplification of " + in[i].get_name());
             simplified_in_self[i] = simplified_curve;
         }
@@ -327,7 +327,7 @@ Clustering_Result one_median_exhaustive(const curve_size_t ell, const Curves &in
 
         for (curve_number_t i = 0; i < in.size(); ++i) {
             Simplification::Subcurve_Shortcut_Graph graph(const_cast<Curve&>(in[i]));
-            auto simplified_curve = graph.weak_minimum_error_simplification(ell);
+            auto simplified_curve = graph.minimum_error_simplification(ell);
             simplified_curve.set_name("Simplification of " + in[i].get_name());
             simplified_in_self[i] = simplified_curve;
         }

include/frechet.hpp

@@ -18,7 +18,7 @@ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLI
 namespace Frechet {
 namespace Continuous {
 
-    extern distance_t epsilon;
+    extern distance_t error;
     extern bool round;
 
     struct Distance {
@@ -39,6 +39,7 @@ namespace Continuous {
             std::vector<std::vector<Interval>>&, std::vector<std::vector<Interval>>&);
 
     distance_t _greedy_upper_bound(const Curve&, const Curve&);
+    distance_t _projective_lower_bound(const Curve&, const Curve&);
 }
 namespace Discrete {
 

include/point.hpp

@@ -15,6 +15,7 @@ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLI
 #include <cmath>
 #include <string>
 #include <sstream>
+#include <type_traits>
 
 #include <pybind11/pybind11.h>
 #include <pybind11/numpy.h>
@@ -92,7 +93,8 @@ class Point : public Coordinates {
         return result;
     }
 
-    inline Point operator*(const distance_t mult) const {
+    template<typename T>
+    inline Point operator*(const T mult) const {
         Point result = *this;
         #pragma omp simd
         for (dimensions_t i = 0; i < dimensions(); ++i){
@@ -146,7 +148,20 @@ class Point : public Coordinates {
         return std::sqrt(length_sqr());
     }
 
-    inline Interval intersection_interval(const distance_t distance_sqr, const Point &line_start, const Point &line_end) const {
+    inline distance_t line_segment_dist_sqr(const Point &p1, const Point &p2) const {
+        const Vector u = p2 - p1;
+        parameter_t projection_param = (*this - p1) * u / (u * u);
+        if (projection_param < parameter_t(0)) projection_param = parameter_t(0);
+        else if (projection_param > parameter_t(1)) projection_param = parameter_t(1);
+        const Point projection = p1 + u * projection_param;
+        return projection.dist_sqr(*this);
+    }
+    
+    inline distance_t line_segment_dist(const Point &p1, const Point &p2) const {
+        return std::sqrt(line_segment_dist(p1, p2));
+    }
+    
+    inline Interval ball_intersection_interval(const distance_t distance_sqr, const Point &line_start, const Point &line_end) const {
         const Vector u = line_end-line_start, v = *this - line_start;
         const parameter_t ulen_sqr = u.length_sqr(), vlen_sqr = v.length_sqr();
 

include/simplification.hpp

@@ -52,7 +52,7 @@ class Subcurve_Shortcut_Graph {
         }
     }
 
-    Curve weak_minimum_error_simplification(const curve_size_t ll) const {
+    Curve minimum_error_simplification(const curve_size_t ll) const {
         if (ll >= curve.complexity()) return curve;
 
         curve_size_t l = ll - 1;
@@ -111,7 +111,7 @@ class Subcurve_Shortcut_Graph {
 
 };
 
-Curve approximate_weak_minimum_link_simplification(const Curve&, const distance_t);
-Curve approximate_weak_minimum_error_simplification(const Curve&, const curve_size_t);
+Curve approximate_minimum_link_simplification(const Curve&, const distance_t);
+Curve approximate_minimum_error_simplification(const Curve&, const curve_size_t);
 
 };

setup.py

@@ -74,7 +74,7 @@ def build_extension(self, ext):
 
 setup(
     name='Fred-Frechet',
-    version='1.7.7',
+    version='1.8',
     author='Dennis Rohde',
     author_email='dennis.rohde@tu-dortmund.de',
     description='A fast, scalable and light-weight C++ Fréchet distance library, exposed to python and focused on (k,l)-clustering of polygonal curves.',

src/curve.cpp

@@ -61,12 +61,12 @@ Curves Curves::simplify(const curve_size_t l, const bool approx = false) {
     Curves result(size(), l, Curves::dimensions());
     for (curve_number_t i = 0; i < size(); ++i) {
         if (approx) {
-            Curve simplified_curve = Simplification::approximate_weak_minimum_error_simplification(std::vector<Curve>::operator[](i), l);
+            Curve simplified_curve = Simplification::approximate_minimum_error_simplification(std::vector<Curve>::operator[](i), l);
             simplified_curve.set_name("Simplification of " + std::vector<Curve>::operator[](i).get_name());
             result[i] = simplified_curve;
         } else {
             Simplification::Subcurve_Shortcut_Graph graph(std::vector<Curve>::operator[](i));
-            Curve simplified_curve = graph.weak_minimum_error_simplification(l);
+            Curve simplified_curve = graph.minimum_error_simplification(l);
             simplified_curve.set_name("Simplification of " + std::vector<Curve>::operator[](i).get_name());
             result[i] = simplified_curve;
         }

src/frechet.cpp

@@ -18,7 +18,7 @@ namespace Frechet {
 
 namespace Continuous {
 
-distance_t epsilon = 0.001;
+distance_t error = 1;
 bool round = true;
 
 std::string Distance::repr() const {
@@ -42,7 +42,7 @@ Distance distance(const Curve &curve1, const Curve &curve2) {
     }
 
     auto start = std::chrono::high_resolution_clock::now();
-    const distance_t lb = std::sqrt(std::max(curve1[0].dist_sqr(curve2[0]), curve1[curve1.complexity()-1].dist_sqr(curve2[curve2.complexity()-1])));
+    const distance_t lb = _projective_lower_bound(curve1, curve2);
     const distance_t ub = _greedy_upper_bound(curve1, curve2);
     auto end = std::chrono::high_resolution_clock::now();
 
@@ -52,7 +52,6 @@ Distance distance(const Curve &curve1, const Curve &curve2) {
 
     auto dist = _distance(curve1, curve2, ub, lb);
     dist.time_bounds = std::chrono::duration_cast<std::chrono::seconds>(end - start).count();
-    if (round) dist.value =  std::round(dist.value * 1e3) / 1e3;
 
     return dist;
 }
@@ -62,9 +61,10 @@ Distance _distance(const Curve &curve1, const Curve &curve2, distance_t ub, dist
     auto start = std::chrono::high_resolution_clock::now();
 
     distance_t split = (ub + lb)/2;
+    const distance_t p_error = lb > 0 ? lb * error / 100 : std::numeric_limits<distance_t>::epsilon();
     std::size_t number_searches = 0;
 
-    if (ub - lb > epsilon) {
+    if (ub - lb > p_error) {
         auto infty = std::numeric_limits<parameter_t>::infinity();
         std::vector<std::vector<parameter_t>> reachable1(curve1.complexity() - 1, std::vector<parameter_t>(curve2.complexity(), infty));
         std::vector<std::vector<parameter_t>> reachable2(curve1.complexity(), std::vector<parameter_t>(curve2.complexity() - 1, infty));
@@ -78,7 +78,7 @@ Distance _distance(const Curve &curve1, const Curve &curve2, distance_t ub, dist
         }
 
         //Binary search over the feasible distances
-        while (ub - lb > epsilon) {
+        while (ub - lb > p_error) {
             ++number_searches;
             split = (ub + lb)/2;
             auto isLessThan = _less_than_or_equal(split, curve1, curve2, reachable1, reachable2, free_intervals1, free_intervals2);
@@ -95,6 +95,10 @@ Distance _distance(const Curve &curve1, const Curve &curve2, distance_t ub, dist
     }
 
     distance_t value = (ub + lb)/2.;
+    const distance_t p_error_exp = std::log10(p_error);
+    const distance_t round_exp = p_error_exp >= 0 ? 0 : static_cast<unsigned int>(-p_error_exp);
+    const distance_t round_fact = std::pow(10, round_exp);
+    if (round) value =  std::round(value * round_fact) / round_fact;
     auto end = std::chrono::high_resolution_clock::now();
     result.value = value;
     result.time_searches = std::chrono::duration_cast<std::chrono::seconds>(end - start).count();
@@ -136,10 +140,10 @@ bool _less_than_or_equal(const distance_t distance, Curve const& curve1, Curve c
     for (curve_size_t i = 0; i < n1; ++i) {
         for (curve_size_t j = 0; j < n2; ++j) {
             if ((i < n1 - 1) and (j > 0)) {
-                free_intervals1[j][i] = curve2[j].intersection_interval(dist_sqr, curve1[i], curve1[i+1]);
+                free_intervals1[j][i] = curve2[j].ball_intersection_interval(dist_sqr, curve1[i], curve1[i+1]);
             }
             if ((j < n2 - 1) and (i > 0)) {
-                free_intervals2[i][j] = curve1[i].intersection_interval(dist_sqr, curve2[j], curve2[j+1]);
+                free_intervals2[i][j] = curve1[i].ball_intersection_interval(dist_sqr, curve2[j], curve2[j+1]);
             }
         }
     }
@@ -201,6 +205,40 @@ distance_t _greedy_upper_bound(const Curve &curve1, const Curve &curve2) {
     return std::sqrt(result);
 }
 
+distance_t _projective_lower_bound(const Curve &curve1, const Curve &curve2) {
+    std::vector<distance_t> distances1_sqr = std::vector<distance_t>(curve2.complexity() - 1), distances2_sqr = std::vector<distance_t>(curve1.complexity() + curve2.complexity() + 2);
+    
+    for (curve_size_t i = 0; i < curve1.complexity(); ++i) {
+        #pragma omp parallel for
+        for (curve_size_t j = 0; j < curve2.complexity() - 1; ++j) {
+            if (curve2[j].dist_sqr(curve2[j+1]) > 0) {
+                distances1_sqr[j] = curve1[i].line_segment_dist_sqr(curve2[j], curve2[j+1]);
+            } else {
+                distances1_sqr[j] = curve1[i].dist_sqr(curve2[j]);
+            }
+        }
+        distances2_sqr[i] = *std::min_element(distances1_sqr.begin(), distances1_sqr.end());
+    }
+    
+    distances1_sqr = std::vector<distance_t>(curve1.complexity() - 1);
+    
+    for (curve_size_t i = 0; i < curve2.complexity(); ++i) {
+        #pragma omp parallel for
+        for (curve_size_t j = 0; j < curve1.complexity() - 1; ++j) {
+            if (curve1[j].dist_sqr(curve1[j+1]) > 0) {
+                distances1_sqr[j] = curve2[i].line_segment_dist_sqr(curve1[j], curve1[j+1]);
+            } else {
+                distances1_sqr[j] = curve2[i].dist_sqr(curve1[j]);
+            }
+        }
+        distances2_sqr[curve1.complexity() + i] = *std::min_element(distances1_sqr.begin(), distances1_sqr.end());
+    }
+    
+    distances2_sqr[curve1.complexity() + curve2.complexity()] = curve1[0].dist_sqr(curve2[0]);
+    distances2_sqr[curve1.complexity() + curve2.complexity() + 1] = curve1[curve1.complexity()-1].dist_sqr(curve2[curve2.complexity()-1]);
+    return std::sqrt(*std::max_element(distances2_sqr.begin(), distances2_sqr.end()));
+}
+
 } // end namespace Continuous
 
 namespace Discrete {