diff --git a/DESCRIPTION b/DESCRIPTION index 677b47b..f295f95 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: mapscanner Title: Print Maps, Draw on Them, Scan Them Back in -Version: 0.1.0.001 +Version: 0.1.0.002 Authors@R: c( person("Mark", "Padgham", , "mark.padgham@email.com", role = c("aut", "cre")), person("Michael D.", "Sumner", , "mdsumner@gmail.com", role = "aut"), diff --git a/codemeta.json b/codemeta.json index 1826f72..3856313 100644 --- a/codemeta.json +++ b/codemeta.json @@ -7,7 +7,7 @@ "codeRepository": "https://github.com/ropensci/mapscanner", "issueTracker": "https://github.com/ropensci/mapscanner/issues", "license": "https://spdx.org/licenses/GPL-3.0", - "version": "0.1.0.001", + "version": "0.1.0.002", "programmingLanguage": { "@type": "ComputerLanguage", "name": "R", diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 9ae0a44..a6ff083 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -33,13 +33,13 @@ BEGIN_RCPP END_RCPP } // rcpp_concaveman -Rcpp::DataFrame rcpp_concaveman(Rcpp::DataFrame xy, Rcpp::IntegerVector hull_in, const double concavity, const double length_threshold); +Rcpp::DataFrame rcpp_concaveman(Rcpp::DataFrame xy, Rcpp::NumericVector hull_in, const double concavity, const double length_threshold); RcppExport SEXP _mapscanner_rcpp_concaveman(SEXP xySEXP, SEXP hull_inSEXP, SEXP concavitySEXP, SEXP length_thresholdSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< Rcpp::DataFrame >::type xy(xySEXP); - Rcpp::traits::input_parameter< Rcpp::IntegerVector >::type hull_in(hull_inSEXP); + Rcpp::traits::input_parameter< Rcpp::NumericVector >::type hull_in(hull_inSEXP); Rcpp::traits::input_parameter< const double >::type concavity(concavitySEXP); Rcpp::traits::input_parameter< const double >::type length_threshold(length_thresholdSEXP); rcpp_result_gen = Rcpp::wrap(rcpp_concaveman(xy, hull_in, concavity, length_threshold)); diff --git a/src/concaveman.cpp b/src/concaveman.cpp index 36f8cf2..d253f7f 100644 --- a/src/concaveman.cpp +++ b/src/concaveman.cpp @@ -3,7 +3,7 @@ //' rcpp_concaveman //' @noRd // [[Rcpp::export]] -Rcpp::DataFrame rcpp_concaveman (Rcpp::DataFrame xy, Rcpp::IntegerVector hull_in, +Rcpp::DataFrame rcpp_concaveman (Rcpp::DataFrame xy, Rcpp::NumericVector hull_in, const double concavity, const double length_threshold) { std::vector x = xy ["x"], y = xy ["y"]; @@ -22,7 +22,7 @@ Rcpp::DataFrame rcpp_concaveman (Rcpp::DataFrame xy, Rcpp::IntegerVector hull_in auto concave_points = concaveman (points, hull, concavity, length_threshold); - Rcpp::IntegerVector xout (concave_points.size ()), + Rcpp::NumericVector xout (concave_points.size ()), yout (concave_points.size ()); for (int i = 0; i < concave_points.size (); i++) { @@ -36,4 +36,3 @@ Rcpp::DataFrame rcpp_concaveman (Rcpp::DataFrame xy, Rcpp::IntegerVector hull_in return res; } - diff --git a/src/concaveman.h b/src/concaveman.h index a5647be..ad15788 100644 --- a/src/concaveman.h +++ b/src/concaveman.h @@ -1,5 +1,8 @@ // // Author: Stanislaw Adaszewski, 2019 +// C++ port from https://github.com/mapbox/concaveman (js) +// +// Comments from js repo added by wheeled // #pragma once @@ -11,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -19,8 +23,11 @@ #include "Rcpp.h" +//#define DEBUG // uncomment to dump debug info to screen +//#define DEBUG_2 // uncomment to dump second-level debug info to screen + template -std::unique_ptr make_unique(Args&&... args) { +std::unique_ptr make_unq_ptr(Args&&... args) { return std::unique_ptr(new T(std::forward(args)...)); } @@ -45,6 +52,7 @@ template T orient2d( } +// check if the edges (p1,q1) and (p2,q2) intersect template bool intersects( const std::array &p1, const std::array &q1, @@ -60,6 +68,7 @@ template bool intersects( } +// square distance between 2 points template T getSqDist( const std::array &p1, const std::array &p2) { @@ -70,6 +79,7 @@ template T getSqDist( } +// square distance from a point to a segment template T sqSegDist( const std::array &p, const std::array &p1, @@ -98,6 +108,7 @@ template T sqSegDist( } +// segment to segment distance, ported from http://geomalgorithms.com/a07-_distance.html by Dan Sunday template T sqSegSegDist(T x0, T y0, T x1, T y1, T x2, T y2, @@ -191,16 +202,16 @@ template class rtree { m_is_leaf(true), m_data(data), m_bounds(bounds) { for (auto i = 0; i < DIM; i++) if (bounds[i] > bounds[i + DIM]) - Rcpp::stop("Bounds minima have to be less than maxima"); // # nocov + throw std::runtime_error("Bounds minima have to be less than maxima"); } void insert(data_type data, const bounds_type &bounds) { if (m_is_leaf) - Rcpp::stop("Cannot insert into leaves"); // # nocov + throw std::runtime_error("Cannot insert into leaves"); m_bounds = updated_bounds(bounds); if (m_children.size() < MAX_CHILDREN) { - auto r = make_unique(data, bounds); + auto r = make_unq_ptr(data, bounds); m_children.push_back(std::move(r)); return; } @@ -216,10 +227,13 @@ template class rtree { } if (!best_child.get().is_leaf()) { best_child.get().insert(data, bounds); +#ifdef DEBUG + std::cout << "best_child: " << bounds[0] << " " << bounds[1] << std::endl; +#endif return; } - auto leaf = make_unique(best_child.get().data(), + auto leaf = make_unq_ptr(best_child.get().data(), best_child.get().bounds()); best_child.get().m_is_leaf = false; best_child.get().m_data = data_type(); @@ -257,7 +271,7 @@ template class rtree { void erase(data_type data, const bounds_type &bounds) { if (m_is_leaf) - Rcpp::stop ("Cannot erase from leaves"); // # nocov + throw std::runtime_error("Cannot erase from leaves"); if (!intersects(bounds)) return; @@ -274,6 +288,23 @@ template class rtree { } } + void print(int level = 0) { + // print the entire tree + + for (auto it = m_children.begin(); it != m_children.end(); ) { + auto bounds = (*it)->m_bounds; + std::string pad(level, '\t'); + if ((*it)->m_is_leaf) { + printf ("%s leaf %0.6f %0.6f \n", pad.c_str(), bounds[0], bounds[1]); + } + else { + printf ("%s branch %0.6f %0.6f %0.6f %0.6f \n", pad.c_str(), bounds[0], bounds[1], bounds[2], bounds[3]); + (*it)->print(level + 1); + } + it++; + } + } + bounds_type updated_bounds(const bounds_type &child_bounds) const { bounds_type res; for (auto i = 0; i < DIM; i++) { @@ -428,8 +459,14 @@ template class CircularList { } ~CircularList() { +#ifdef DEBUG + std::cout << "~CircularList()" << std::endl; +#endif auto node = m_last; while (true) { +#ifdef DEBUG +// std::cout << (i++) << std::endl; +#endif auto tmp = node; node = node->m_next; delete tmp; @@ -442,7 +479,7 @@ template class CircularList { auto elem = new CircularElement(std::forward(args)...); if (prev == nullptr && m_last != nullptr) - Rcpp::stop ("Once the list is non-empty you must specify where to insert"); // # nocov + throw std::runtime_error("Once the list is non-empty you must specify where to insert"); if (prev == nullptr) { elem->m_prev = elem->m_next = elem; @@ -464,6 +501,7 @@ template class CircularList { }; +// update the bounding box of a node's edge template void updateBBox(typename CircularElement::ptr_type elem) { auto &node(elem->data()); auto p1 = node.p; @@ -475,10 +513,40 @@ template void updateBBox(typename CircularElement::ptr_type elem) { } +#ifdef DEBUG_2 +template void snapshot( + const std::array &a, + const std::array &b, + const std::array &c, + const std::array &d, + const double sqLen, + const double maxSqLen, + const std::array &trigger, + const bool use_trigger) { + + if ( !use_trigger || trigger == b ) { + if ( !use_trigger ) + printf ("Snapshot untriggered\n"); + else + printf ("Snapshot trigger: %0.6f %0.6f \n", trigger[0], trigger[1]); + printf ("... segment a, b: %0.6f %0.6f, %0.6f %0.6f \n", a[0], a[1], b[0], b[1]); + printf ("... segment c, d: %0.6f %0.6f, %0.6f %0.6f \n", c[0], c[1], d[0], d[1]); + printf ("... sqDist a-b, b-c, c-d: %e, %e, %e", getSqDist(a, b), getSqDist(b, c), getSqDist(c, d)); + printf ("... sqLen, maxSqLen: %e, %e", sqLen, maxSqLen); + } +} +#endif + + template std::vector> concaveman( const std::vector> &points, + // start with a convex hull of the points const std::vector &hull, - T concavity=2, T lengthThreshold=0) { + // a relative measure of concavity; higher value means simpler hull + T concavity=2, + // when a segment goes below this length threshold, it won't be drilled down further + T lengthThreshold=0 + ) { typedef Node node_type; typedef std::array point_type; @@ -486,13 +554,18 @@ template std::vector> concaveman( typedef CircularList circ_list_type; typedef circ_elem_type *circ_elem_ptr_type; +#ifdef DEBUG + std::cout << "concaveman()" << std::endl; +#endif + // exit if hull includes all points already if (hull.size() == points.size()) { - std::vector res; // # nocov - for (auto &i : hull) res.push_back(points[i]); // # nocov - return res; // # nocov + std::vector res; + for (auto &i : hull) res.push_back(points[i]); + return res; } + // index the points with an R-tree rtree tree; for (auto &p : points) @@ -503,6 +576,7 @@ template std::vector> concaveman( std::list queue; + // turn the convex hull into a linked list and populate the initial edge queue with the nodes for (auto &idx : hull) { auto &p = points[idx]; tree.erase(p, { p[0], p[1], p[0], p[1] }); @@ -510,11 +584,23 @@ template std::vector> concaveman( queue.push_back(last); } +#ifdef DEBUG_2 + tree.print(0); +#endif + + // loops through the hull? why? +#ifdef DEBUG + std::cout << "Starting hull: "; +#endif for (auto elem = last->next(); ; elem=elem->next()) { +#ifdef DEBUG + std::cout << elem->data().p[0] << " " << elem->data().p[1] << std::endl; +#endif if (elem == last) break; } + // index the segments with an R-tree (for intersection checks) rtree segTree; for (auto &elem : queue) { auto &node(elem->data()); @@ -526,28 +612,45 @@ template std::vector> concaveman( auto sqConcavity = concavity * concavity; auto sqLenThreshold = lengthThreshold * lengthThreshold; + // process edges one by one while (!queue.empty()) { auto elem = *queue.begin(); queue.pop_front(); - auto a = elem->data().p; - auto b = elem->next()->data().p; + auto a = elem->prev()->data().p; + auto b = elem->data().p; + auto c = elem->next()->data().p; + auto d = elem->next()->next()->data().p; - auto sqLen = getSqDist(a, b); + // skip the edge if it's already short enough + auto sqLen = getSqDist(b, c); if (sqLen < sqLenThreshold) continue; auto maxSqLen = sqLen / sqConcavity; +#ifdef DEBUG_2 + // dump key parameters either on every pass or when a certain point is 'b' + point_type trigger = { 151.1373474787800, -33.7733192376544 }; + snapshot(a, b, c, d, sqLen, maxSqLen, trigger, true); +#endif + + // find the best connection point for the current edge to flex inward to bool ok; - auto p = findCandidate(tree, elem->prev()->data().p, a, b, - elem->next()->next()->data().p, maxSqLen, segTree, ok); + auto p = findCandidate(tree, a, b, c, d, maxSqLen, segTree, ok); + + // if we found a connection and it satisfies our concavity measure + if (ok && std::min(getSqDist(p, b), getSqDist(p, c)) <= maxSqLen) { - if (ok && std::min(getSqDist(p, a), getSqDist(p, b)) <= maxSqLen) { +#ifdef DEBUG + printf ("Modifying hull, p: %0.6f %0.6f \n" ,p[0], p[1]); +#endif + // connect the edge endpoints through this point and add 2 new edges to the queue queue.push_back(elem); queue.push_back(elem->insert(p)); + // update point and segment indexes auto &node = elem->data(); auto &next = elem->next()->data(); @@ -560,8 +663,13 @@ template std::vector> concaveman( segTree.insert(elem, { node.minX, node.minY, node.maxX, node.maxY }); segTree.insert(elem->next(), { next.minX, next.minY, next.maxX, next.maxY }); } +#ifdef DEBUG + else + printf ("No point found along segment: %0.6f %0.6f, %0.6f %0.6f \n", b[0], b[1], c[0], c[1]); +#endif } + // convert the resulting hull linked list to an array of points std::vector concave; for (auto elem = last->next(); ; elem = elem->next()) { concave.push_back(elem->data().p); @@ -584,17 +692,23 @@ template std::array findCandidate( bool &ok) { typedef std::array point_type; - //typedef CircularElement> circ_elem_type; + typedef CircularElement> circ_elem_type; typedef rtree> tree_type; typedef const tree_type const_tree_type; typedef std::reference_wrapper tree_ref_type; typedef std::tuple tuple_type; +#ifdef DEBUG + std::cout << "findCandidate(), maxDist: " << maxDist << std::endl; +#endif + ok = false; std::priority_queue, compare_first> queue; std::reference_wrapper node = tree; + // search through the point R-tree with a depth-first search using a priority queue + // in the order of distance to the edge (b, c) while (true) { for (auto &child : node.get().children()) { @@ -603,7 +717,7 @@ template std::array findCandidate( auto dist = child->is_leaf() ? sqSegDist(pt, b, c) : sqSegBoxDist(b, c, *child); if (dist > maxDist) - continue; + continue; // skip the node if it's farther than we ever need queue.push(tuple_type(-dist, *child)); } @@ -615,29 +729,46 @@ template std::array findCandidate( auto bounds = std::get<1>(item).get().bounds(); point_type p = { bounds[0], bounds[1] }; + // skip all points that are as close to adjacent edges (a,b) and (c,d), + // and points that would introduce self-intersections when connected auto d0 = sqSegDist(p, a, b); auto d1 = sqSegDist(p, c, d); - if (std::get<0>(item) < d0 && std::get<0>(item) < d1 && +#ifdef DEBUG_2 + printf (" p: %0.6f %0.6f sqSegDist: %e, %e, %e \n", bounds[0], bounds[1], d0, std::get<0>(item), d1); +#endif + + if (-std::get<0>(item) < d0 && -std::get<0>(item) < d1 && noIntersections(b, p, segTree) && noIntersections(c, p, segTree)) { ok = true; return std::get<1>(item).get().data(); } - } + +#ifdef DEBUG_2 + else { + bool cond1 = -std::get<0>(item) < d0; + bool cond2 = -std::get<0>(item) < d1; + bool cond3 = noIntersections(b, p, segTree); + bool cond4 = noIntersections(c, p, segTree); + std::cout << "Not OK: " << cond1 << " " << cond2 << " " << cond3 << " " << cond4 << std::endl; + } +#endif + } if (queue.empty()) - break; // # nocov + break; node = std::get<1>(queue.top()); queue.pop(); } - return point_type(); // # nocov + return point_type(); } +// square distance from a segment bounding box to the given one template T sqSegBoxDist( const std::array &a, const std::array &b, @@ -687,6 +818,7 @@ template bool inside( } +// check if the edge (a,b) doesn't intersect any other edges template bool noIntersections( const std::array &a, const std::array &b, @@ -703,11 +835,8 @@ template bool noIntersections( auto elem = ch.data(); if (intersects(elem->data().p, elem->next()->data().p, a, b)) - return false; // # nocov + return false; } return true; } - -Rcpp::DataFrame rcpp_concaveman (Rcpp::DataFrame xy, Rcpp::IntegerVector hull_in, - const double concavity, const double length_threshold);