Feature: AABB tree for lines

Small refactoring of AABB tree distance query function, to allow different primitives (apart from triangles)
Implemented Distancer and functions to create AABB tree from lines and use closest point query
Added test for the AABBTree with lines
Added Benchmark comparing EdgeGrid with AABBTree on line contours (Inside AABBTree test file, disabled under compilation flag)

src/libslic3r/AABBTreeIndirect.hpp

@@ -446,21 +446,6 @@ namespace detail {
 		}
 	}
 
-	// Nothing to do with COVID-19 social distancing.
-	template<typename AVertexType, typename AIndexedFaceType, typename ATreeType, typename AVectorType>
-	struct IndexedTriangleSetDistancer {
-		using VertexType 		= AVertexType;
-		using IndexedFaceType 	= AIndexedFaceType;
-		using TreeType			= ATreeType;
-		using VectorType 		= AVectorType;
-
-		const std::vector<VertexType> 		&vertices;
-		const std::vector<IndexedFaceType> 	&faces;
-		const TreeType 						&tree;
-
-		const VectorType					 origin;
-	};
-
     // Real-time collision detection, Ericson, Chapter 5
     template<typename Vector>
     static inline Vector closest_point_to_triangle(const Vector &p, const Vector &a, const Vector &b, const Vector &c)
@@ -511,16 +496,44 @@ namespace detail {
         return a + ab * v + ac * w; // = u*a + v*b + w*c, u = va * denom = 1.0-v-w
     };
 
-	template<typename IndexedTriangleSetDistancerType, typename Scalar>
+
-    static inline Scalar squared_distance_to_indexed_triangle_set_recursive(
+	// Nothing to do with COVID-19 social distancing.
-        IndexedTriangleSetDistancerType	&distancer,
+	template<typename AVertexType, typename AIndexedFaceType, typename ATreeType, typename AVectorType>
+	struct IndexedTriangleSetDistancer {
+		using VertexType 		= AVertexType;
+		using IndexedFaceType 	= AIndexedFaceType;
+		using TreeType			= ATreeType;
+		using VectorType 		= AVectorType;
+		using ScalarType 		= typename VectorType::Scalar;
+
+		const std::vector<VertexType> 		&vertices;
+		const std::vector<IndexedFaceType> 	&faces;
+		const TreeType 						&tree;
+
+		const VectorType					 origin;
+
+		inline VectorType closest_point_to_origin(size_t primitive_index,
+		        ScalarType& squared_distance){
+		    const auto &triangle = this->faces[primitive_index];
+		    VectorType closest_point = closest_point_to_triangle<VectorType>(origin,
+		            this->vertices[triangle(0)].template cast<ScalarType>(),
+		            this->vertices[triangle(1)].template cast<ScalarType>(),
+		            this->vertices[triangle(2)].template cast<ScalarType>());
+		    squared_distance = (origin - closest_point).squaredNorm();
+		    return closest_point;
+		}
+	};
+
+	template<typename IndexedPrimitivesDistancerType, typename Scalar>
+    static inline Scalar squared_distance_to_indexed_primitives_recursive(
+        IndexedPrimitivesDistancerType	&distancer,
 		size_t 							 node_idx,
 		Scalar 							 low_sqr_d,
   		Scalar 							 up_sqr_d,
 		size_t 							&i,
-  		Eigen::PlainObjectBase<typename IndexedTriangleSetDistancerType::VectorType> &c)
+  		Eigen::PlainObjectBase<typename IndexedPrimitivesDistancerType::VectorType> &c)
 	{
-		using Vector = typename IndexedTriangleSetDistancerType::VectorType;
+		using Vector = typename IndexedPrimitivesDistancerType::VectorType;
 
   		if (low_sqr_d > up_sqr_d)
 			return low_sqr_d;
@@ -538,13 +551,9 @@ namespace detail {
 		assert(node.is_valid());
   		if (node.is_leaf()) 
   		{
-            const auto &triangle = distancer.faces[node.idx];
+            Scalar sqr_dist;
-            Vector c_candidate = closest_point_to_triangle<Vector>(
+            Vector c_candidate = distancer.closest_point_to_origin(node.idx, sqr_dist);
-				distancer.origin, 
+            set_min(sqr_dist, node.idx, c_candidate);
-                distancer.vertices[triangle(0)].template cast<Scalar>(),
-                distancer.vertices[triangle(1)].template cast<Scalar>(),
-                distancer.vertices[triangle(2)].template cast<Scalar>());
-            set_min((c_candidate - distancer.origin).squaredNorm(), node.idx, c_candidate);
   		} 
   		else
   		{
@@ -561,7 +570,7 @@ namespace detail {
 			{
                 size_t	i_left;
                 Vector 	c_left = c;
-                Scalar	sqr_d_left = squared_distance_to_indexed_triangle_set_recursive(distancer, left_node_idx, low_sqr_d, up_sqr_d, i_left, c_left);
+                Scalar	sqr_d_left = squared_distance_to_indexed_primitives_recursive(distancer, left_node_idx, low_sqr_d, up_sqr_d, i_left, c_left);
 				set_min(sqr_d_left, i_left, c_left);
 				looked_left = true;
 			};
@@ -569,13 +578,13 @@ namespace detail {
 			{
                 size_t	i_right;
                 Vector	c_right = c;
-                Scalar	sqr_d_right = squared_distance_to_indexed_triangle_set_recursive(distancer, right_node_idx, low_sqr_d, up_sqr_d, i_right, c_right);
+                Scalar	sqr_d_right = squared_distance_to_indexed_primitives_recursive(distancer, right_node_idx, low_sqr_d, up_sqr_d, i_right, c_right);
 				set_min(sqr_d_right, i_right, c_right);
 				looked_right = true;
 			};
 
 			// must look left or right if in box
-            using BBoxScalar = typename IndexedTriangleSetDistancerType::TreeType::BoundingBox::Scalar;
+            using BBoxScalar = typename IndexedPrimitivesDistancerType::TreeType::BoundingBox::Scalar;
             if (node_left.bbox.contains(distancer.origin.template cast<BBoxScalar>()))
 			  	look_left();
             if (node_right.bbox.contains(distancer.origin.template cast<BBoxScalar>()))
@@ -747,7 +756,7 @@ inline typename VectorType::Scalar squared_distance_to_indexed_triangle_set(
     auto distancer = detail::IndexedTriangleSetDistancer<VertexType, IndexedFaceType, TreeType, VectorType>
         { vertices, faces, tree, point };
     return tree.empty() ? Scalar(-1) : 
-    	detail::squared_distance_to_indexed_triangle_set_recursive(distancer, size_t(0), Scalar(0), std::numeric_limits<Scalar>::infinity(), hit_idx_out, hit_point_out);
+    	detail::squared_distance_to_indexed_primitives_recursive(distancer, size_t(0), Scalar(0), std::numeric_limits<Scalar>::infinity(), hit_idx_out, hit_point_out);
 }
 
 // Decides if exists some triangle in defined radius on a 3D indexed triangle set using a pre-built AABBTreeIndirect::Tree.
@@ -779,7 +788,7 @@ inline bool is_any_triangle_in_radius(
 		return false;
 	}
 
-	detail::squared_distance_to_indexed_triangle_set_recursive(distancer, size_t(0), Scalar(0), max_distance_squared, hit_idx, hit_point);
+	detail::squared_distance_to_indexed_primitives_recursive(distancer, size_t(0), Scalar(0), max_distance_squared, hit_idx, hit_point);
 
     return hit_point.allFinite();
 }

src/libslic3r/AABBTreeLines.hpp

@@ -0,0 +1,112 @@
+#ifndef SRC_LIBSLIC3R_AABBTREELINES_HPP_
+#define SRC_LIBSLIC3R_AABBTREELINES_HPP_
+
+#include "libslic3r/Point.hpp"
+#include "libslic3r/EdgeGrid.hpp"
+#include "libslic3r/AABBTreeIndirect.hpp"
+#include "libslic3r/Line.hpp"
+
+namespace Slic3r {
+
+namespace AABBTreeLines {
+
+namespace detail {
+
+template<typename ALineType, typename ATreeType, typename AVectorType>
+struct IndexedLinesDistancer {
+    using LineType = ALineType;
+    using TreeType = ATreeType;
+    using VectorType = AVectorType;
+    using ScalarType = typename VectorType::Scalar;
+
+    const std::vector<LineType> &lines;
+    const TreeType &tree;
+
+    const VectorType origin;
+
+    inline VectorType closest_point_to_origin(size_t primitive_index,
+            ScalarType &squared_distance) {
+        VectorType nearest_point;
+        const LineType &line = lines[primitive_index];
+        squared_distance = line_alg::distance_to_squared(line, origin, &nearest_point);
+        return nearest_point;
+    }
+};
+
+}
+
+// Build a balanced AABB Tree over a vector of float lines, balancing the tree
+// on centroids of the lines.
+// Epsilon is applied to the bounding boxes of the AABB Tree to cope with numeric inaccuracies
+// during tree traversal.
+template<typename LineType>
+inline AABBTreeIndirect::Tree<2, typename LineType::Scalar> build_aabb_tree_over_indexed_lines(
+        const std::vector<LineType> &lines,
+        //FIXME do we want to apply an epsilon?
+        const float eps = 0)
+        {
+    using TreeType = AABBTreeIndirect::Tree<2, typename LineType::Scalar>;
+//    using              CoordType      = typename TreeType::CoordType;
+    using VectorType = typename TreeType::VectorType;
+    using BoundingBox = typename TreeType::BoundingBox;
+
+    struct InputType {
+        size_t idx() const {
+            return m_idx;
+        }
+        const BoundingBox& bbox() const {
+            return m_bbox;
+        }
+        const VectorType& centroid() const {
+            return m_centroid;
+        }
+
+        size_t m_idx;
+        BoundingBox m_bbox;
+        VectorType m_centroid;
+    };
+
+    std::vector<InputType> input;
+    input.reserve(lines.size());
+    const VectorType veps(eps, eps);
+    for (size_t i = 0; i < lines.size(); ++i) {
+        const LineType &line = lines[i];
+        InputType n;
+        n.m_idx = i;
+        n.m_centroid = (line.a + line.b) * 0.5;
+        n.m_bbox = BoundingBox(line.a, line.a);
+        n.m_bbox.extend(line.b);
+        n.m_bbox.min() -= veps;
+        n.m_bbox.max() += veps;
+        input.emplace_back(n);
+    }
+
+    TreeType out;
+    out.build(std::move(input));
+    return out;
+}
+
+// Finding a closest line, its closest point and squared distance to the closest point
+// Returns squared distance to the closest point or -1 if the input is empty.
+template<typename LineType, typename TreeType, typename VectorType>
+inline typename VectorType::Scalar squared_distance_to_indexed_lines(
+        const std::vector<LineType> &lines,
+        const TreeType &tree,
+        const VectorType &point,
+        size_t &hit_idx_out,
+        Eigen::PlainObjectBase<VectorType> &hit_point_out)
+        {
+    using Scalar = typename VectorType::Scalar;
+    auto distancer = detail::IndexedLinesDistancer<LineType, TreeType, VectorType>
+            { lines, tree, point };
+    return tree.empty() ?
+                          Scalar(-1) :
+                          AABBTreeIndirect::detail::squared_distance_to_indexed_primitives_recursive(distancer, size_t(0), Scalar(0),
+                                  std::numeric_limits<Scalar>::infinity(), hit_idx_out, hit_point_out);
+}
+
+}
+
+}
+
+#endif /* SRC_LIBSLIC3R_AABBTREELINES_HPP_ */

tests/libslic3r/test_aabbindirect.cpp

@@ -3,6 +3,7 @@
 
 #include <libslic3r/TriangleMesh.hpp>
 #include <libslic3r/AABBTreeIndirect.hpp>
+#include <libslic3r/AABBTreeLines.hpp>
 
 using namespace Slic3r;
 
@@ -58,3 +59,287 @@ TEST_CASE("Building a tree over a box, ray caster and closest query", "[AABBIndi
     REQUIRE(closest_point.y() == Approx(0.5));
     REQUIRE(closest_point.z() == Approx(1.));
 }
+
+TEST_CASE("Creating a several 2d lines, testing closest point query", "[AABBIndirect]")
+{
+    std::vector<Linef> lines { };
+    lines.push_back(Linef(Vec2d(0.0, 0.0), Vec2d(1.0, 0.0)));
+    lines.push_back(Linef(Vec2d(1.0, 0.0), Vec2d(1.0, 1.0)));
+    lines.push_back(Linef(Vec2d(1.0, 1.0), Vec2d(0.0, 1.0)));
+    lines.push_back(Linef(Vec2d(0.0, 1.0), Vec2d(0.0, 0.0)));
+
+    auto tree = AABBTreeLines::build_aabb_tree_over_indexed_lines(lines);
+
+    size_t hit_idx_out;
+    Vec2d hit_point_out;
+    auto sqr_dist = AABBTreeLines::squared_distance_to_indexed_lines(lines, tree, Vec2d(0.0, 0.0), hit_idx_out,
+            hit_point_out);
+    REQUIRE(sqr_dist == Approx(0.0));
+    REQUIRE((hit_idx_out == 0 || hit_idx_out == 3));
+    REQUIRE(hit_point_out.x() == Approx(0.0));
+    REQUIRE(hit_point_out.y() == Approx(0.0));
+
+    sqr_dist = AABBTreeLines::squared_distance_to_indexed_lines(lines, tree, Vec2d(1.5, 0.5), hit_idx_out,
+            hit_point_out);
+    REQUIRE(sqr_dist == Approx(0.25));
+    REQUIRE(hit_idx_out == 1);
+    REQUIRE(hit_point_out.x() == Approx(1.0));
+    REQUIRE(hit_point_out.y() == Approx(0.5));
+}
+
+#if 0
+#include "libslic3r/EdgeGrid.hpp"
+#include <iostream>
+#include <ctime>
+#include <ratio>
+#include <chrono>
+TEST_CASE("AABBTreeLines vs SignedDistanceGrid time Benchmark", "[AABBIndirect]")
+{
+    std::vector<Points> lines { Points { } };
+    std::vector<Linef> linesf { };
+    Vec2d prevf { };
+
+    // NOTE: max coord value of the lines is approx 83 mm
+    for (int r = 1; r < 1000; ++r) {
+        lines[0].push_back(Point::new_scale(Vec2d(exp(0.005f * r) * cos(r), exp(0.005f * r) * cos(r))));
+        linesf.emplace_back(prevf, Vec2d(exp(0.005f * r) * cos(r), exp(0.005f * r) * cos(r)));
+        prevf = linesf.back().b;
+    }
+
+    int build_num = 10000;
+    using namespace std::chrono;
+    {
+        std::cout << "building the tree " << build_num << " times..." << std::endl;
+        high_resolution_clock::time_point t1 = high_resolution_clock::now();
+        for (int i = 0; i < build_num; ++i) {
+            volatile auto tree = AABBTreeLines::build_aabb_tree_over_indexed_lines(linesf);
+        }
+        high_resolution_clock::time_point t2 = high_resolution_clock::now();
+
+        duration<double> time_span = duration_cast<duration<double>>(t2 - t1);
+        std::cout << "It took " << time_span.count() << " seconds." << std::endl << std::endl;
+
+    }
+    {
+        std::cout << "building the grid res 1mm ONLY " << build_num/100 << " !!! times..." << std::endl;
+        high_resolution_clock::time_point t1 = high_resolution_clock::now();
+        for (int i = 0; i < build_num/100; ++i) {
+            EdgeGrid::Grid grid { };
+            grid.create(lines, scaled(1.0), true);
+            grid.calculate_sdf();
+        }
+        high_resolution_clock::time_point t2 = high_resolution_clock::now();
+        duration<double> time_span = duration_cast<duration<double>>(t2 - t1);
+        std::cout << "It took " << time_span.count() << " seconds." << std::endl << std::endl;
+    }
+    {
+        std::cout << "building the grid res 10mm " << build_num << " times..." << std::endl;
+        high_resolution_clock::time_point t1 = high_resolution_clock::now();
+        for (int i = 0; i < build_num; ++i) {
+            EdgeGrid::Grid grid { };
+            grid.create(lines, scaled(10.0), true);
+            grid.calculate_sdf();
+        }
+        high_resolution_clock::time_point t2 = high_resolution_clock::now();
+        duration<double> time_span = duration_cast<duration<double>>(t2 - t1);
+        std::cout << "It took " << time_span.count() << " seconds." << std::endl << std::endl;
+    }
+
+    EdgeGrid::Grid grid10 { };
+    grid10.create(lines, scaled(10.0), true);
+    coord_t query10_res = scaled(10.0);
+    grid10.calculate_sdf();
+
+    EdgeGrid::Grid grid1 { };
+    grid1.create(lines, scaled(1.0), true);
+    coord_t query1_res = scaled(1.0);
+    grid1.calculate_sdf();
+
+    auto tree = AABBTreeLines::build_aabb_tree_over_indexed_lines(linesf);
+
+    int query_num = 10000;
+    Points query_points { };
+    std::vector<Vec2d> query_pointsf { };
+    for (int x = 0; x < query_num; ++x) {
+        Vec2d qp { rand() / (double(RAND_MAX) + 1.0f) * 200.0 - 100.0, rand() / (double(RAND_MAX) + 1.0f) * 200.0
+                - 100.0 };
+        query_pointsf.push_back(qp);
+        query_points.push_back(Point::new_scale(qp));
+    }
+
+    {
+        std::cout << "querying tree " << query_num << " times..." << std::endl;
+        high_resolution_clock::time_point t1 = high_resolution_clock::now();
+        for (const Vec2d &qp : query_pointsf) {
+            size_t hit_idx_out;
+            Vec2d hit_point_out;
+            AABBTreeLines::squared_distance_to_indexed_lines(linesf, tree, qp, hit_idx_out, hit_point_out);
+        }
+        high_resolution_clock::time_point t2 = high_resolution_clock::now();
+        duration<double> time_span = duration_cast<duration<double>>(t2 - t1);
+        std::cout << "It took " << time_span.count() << " seconds." << std::endl << std::endl;
+    }
+
+    {
+        std::cout << "querying grid res 1mm " << query_num << " times..." << std::endl;
+        high_resolution_clock::time_point t1 = high_resolution_clock::now();
+        for (const Point &qp : query_points) {
+            volatile auto dist = grid1.closest_point_signed_distance(qp, query1_res);
+        }
+        high_resolution_clock::time_point t2 = high_resolution_clock::now();
+        duration<double> time_span = duration_cast<duration<double>>(t2 - t1);
+        std::cout << "It took " << time_span.count() << " seconds." << std::endl << std::endl;
+    }
+
+    {
+        std::cout << "querying grid res 10mm " << query_num << " times..." << std::endl;
+        high_resolution_clock::time_point t1 = high_resolution_clock::now();
+        for (const Point &qp : query_points) {
+            volatile auto dist = grid10.closest_point_signed_distance(qp, query10_res);
+        }
+        high_resolution_clock::time_point t2 = high_resolution_clock::now();
+        duration<double> time_span = duration_cast<duration<double>>(t2 - t1);
+        std::cout << "It took " << time_span.count() << " seconds." << std::endl << std::endl;
+    }
+
+    std::cout << "Test build and queries together - same number of contour points and query points" << std::endl << std::endl;
+
+    std::vector<int> point_counts { 100, 300, 500, 1000, 3000 };
+    for (auto count : point_counts) {
+
+        std::vector<Points> lines { Points { } };
+        std::vector<Linef> linesf { };
+        Vec2d prevf { };
+        Points query_points { };
+        std::vector<Vec2d> query_pointsf { };
+
+        for (int x = 0; x < count; ++x) {
+            Vec2d cp { rand() / (double(RAND_MAX) + 1.0f) * 200.0 - 100.0, rand() / (double(RAND_MAX) + 1.0f) * 200.0
+                    - 100.0 };
+            lines[0].push_back(Point::new_scale(cp));
+            linesf.emplace_back(prevf, cp);
+            prevf = linesf.back().b;
+
+            Vec2d qp { rand() / (double(RAND_MAX) + 1.0f) * 200.0 - 100.0, rand() / (double(RAND_MAX) + 1.0f) * 200.0
+                    - 100.0 };
+            query_pointsf.push_back(qp);
+            query_points.push_back(Point::new_scale(qp));
+        }
+
+        std::cout << "Test for point count: " << count << std::endl;
+        {
+            high_resolution_clock::time_point t1 = high_resolution_clock::now();
+            auto tree = AABBTreeLines::build_aabb_tree_over_indexed_lines(linesf);
+            for (const Vec2d &qp : query_pointsf) {
+                size_t hit_idx_out;
+                Vec2d hit_point_out;
+                AABBTreeLines::squared_distance_to_indexed_lines(linesf, tree, qp, hit_idx_out, hit_point_out);
+            }
+            high_resolution_clock::time_point t2 = high_resolution_clock::now();
+            duration<double> time_span = duration_cast<duration<double>>(t2 - t1);
+            std::cout << "    Tree     took " << time_span.count() << " seconds." << std::endl;
+        }
+
+        {
+            high_resolution_clock::time_point t1 = high_resolution_clock::now();
+            EdgeGrid::Grid grid1 { };
+            grid1.create(lines, scaled(1.0), true);
+            coord_t query1_res = scaled(1.0);
+            grid1.calculate_sdf();
+            for (const Point &qp : query_points) {
+                volatile auto dist = grid1.closest_point_signed_distance(qp, query1_res);
+            }
+            high_resolution_clock::time_point t2 = high_resolution_clock::now();
+            duration<double> time_span = duration_cast<duration<double>>(t2 - t1);
+            std::cout << "    Grid 1mm  took " << time_span.count() << " seconds." << std::endl;
+        }
+
+        {
+            high_resolution_clock::time_point t1 = high_resolution_clock::now();
+            EdgeGrid::Grid grid10 { };
+            grid10.create(lines, scaled(10.0), true);
+            coord_t query10_res = scaled(10.0);
+            grid10.calculate_sdf();
+            for (const Point &qp : query_points) {
+                volatile auto dist = grid10.closest_point_signed_distance(qp, query10_res);
+            }
+            high_resolution_clock::time_point t2 = high_resolution_clock::now();
+            duration<double> time_span = duration_cast<duration<double>>(t2 - t1);
+            std::cout << "    Grid 10mm took " << time_span.count() << " seconds." << std::endl;
+        }
+
+    }
+
+
+
+    std::cout << "Test build and queries together - same number of contour points and query points" << std::endl <<
+            "And with limited contour edge length to 4mm " << std::endl;
+      for (auto count : point_counts) {
+
+          std::vector<Points> lines { Points { } };
+          std::vector<Linef> linesf { };
+          Vec2d prevf { };
+          Points query_points { };
+          std::vector<Vec2d> query_pointsf { };
+
+          for (int x = 0; x < count; ++x) {
+              Vec2d cp { rand() / (double(RAND_MAX) + 1.0f) * 200.0 - 100.0, rand() / (double(RAND_MAX) + 1.0f) * 200.0
+                      - 100.0 };
+              Vec2d contour = prevf + cp.normalized()*4.0; // limits the cnotour edge len to 4mm
+              lines[0].push_back(Point::new_scale(contour));
+              linesf.emplace_back(prevf, contour);
+              prevf = linesf.back().b;
+
+              Vec2d qp { rand() / (double(RAND_MAX) + 1.0f) * 200.0 - 100.0, rand() / (double(RAND_MAX) + 1.0f) * 200.0
+                      - 100.0 };
+              query_pointsf.push_back(qp);
+              query_points.push_back(Point::new_scale(qp));
+          }
+
+          std::cout << "Test for point count: " << count << std::endl;
+          {
+              high_resolution_clock::time_point t1 = high_resolution_clock::now();
+              auto tree = AABBTreeLines::build_aabb_tree_over_indexed_lines(linesf);
+              for (const Vec2d &qp : query_pointsf) {
+                  size_t hit_idx_out;
+                  Vec2d hit_point_out;
+                  AABBTreeLines::squared_distance_to_indexed_lines(linesf, tree, qp, hit_idx_out, hit_point_out);
+              }
+              high_resolution_clock::time_point t2 = high_resolution_clock::now();
+              duration<double> time_span = duration_cast<duration<double>>(t2 - t1);
+              std::cout << "    Tree     took " << time_span.count() << " seconds." << std::endl;
+          }
+
+          {
+              high_resolution_clock::time_point t1 = high_resolution_clock::now();
+              EdgeGrid::Grid grid1 { };
+              grid1.create(lines, scaled(1.0), true);
+              coord_t query1_res = scaled(1.0);
+              grid1.calculate_sdf();
+              for (const Point &qp : query_points) {
+                  volatile auto dist = grid1.closest_point_signed_distance(qp, query1_res);
+              }
+              high_resolution_clock::time_point t2 = high_resolution_clock::now();
+              duration<double> time_span = duration_cast<duration<double>>(t2 - t1);
+              std::cout << "    Grid 1mm  took " << time_span.count() << " seconds." << std::endl;
+          }
+
+          {
+              high_resolution_clock::time_point t1 = high_resolution_clock::now();
+              EdgeGrid::Grid grid10 { };
+              grid10.create(lines, scaled(10.0), true);
+              coord_t query10_res = scaled(10.0);
+              grid10.calculate_sdf();
+              for (const Point &qp : query_points) {
+                  volatile auto dist = grid10.closest_point_signed_distance(qp, query10_res);
+              }
+              high_resolution_clock::time_point t2 = high_resolution_clock::now();
+              duration<double> time_span = duration_cast<duration<double>>(t2 - t1);
+              std::cout << "    Grid 10mm took " << time_span.count() << " seconds." << std::endl;
+          }
+
+      }
+
+}
+#endif
+