From 143d1e57cf0c5f056eb6fe6e7a84aa02a63e8879 Mon Sep 17 00:00:00 2001 From: Artem Pavlenko Date: Thu, 21 Sep 2017 15:40:30 +0200 Subject: [PATCH] add 'closest_point' geometry algorithm --- .../extensions/algorithms/closest_point.hpp | 221 ++++++++++++++++++ .../strategies/cartesian/closest_point.hpp | 164 +++++++++++++ include/mapnik/geometry/closest_point.hpp | 144 ++++++++++++ test/unit/geometry/closest_point.cpp | 46 ++++ 4 files changed, 575 insertions(+) create mode 100644 deps/boost/geometry/extensions/algorithms/closest_point.hpp create mode 100644 deps/boost/geometry/extensions/strategies/cartesian/closest_point.hpp create mode 100644 include/mapnik/geometry/closest_point.hpp create mode 100644 test/unit/geometry/closest_point.cpp diff --git a/deps/boost/geometry/extensions/algorithms/closest_point.hpp b/deps/boost/geometry/extensions/algorithms/closest_point.hpp new file mode 100644 index 000000000..97deb2d57 --- /dev/null +++ b/deps/boost/geometry/extensions/algorithms/closest_point.hpp @@ -0,0 +1,221 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2007-2013 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2008-2013 Bruno Lalande, Paris, France. +// Copyright (c) 2009-2013 Mateusz Loskot, London, UK. +// Copyright (c) 2013 Adam Wulkiewicz, Lodz, Poland. +// Copyright (c) 2017 Artem Pavlenko, UK + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_ALGORITHMS_CLOSEST_POINT_HPP +#define BOOST_GEOMETRY_ALGORITHMS_CLOSEST_POINT_HPP + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include + +#include +#include + +#include +#include + +#include + + +namespace boost { namespace geometry +{ + +#ifndef DOXYGEN_NO_DETAIL +namespace detail { namespace closest_point { + + +template +struct point_point +{ + template + static inline void apply(Point1 const& point1, Point2 const& point2, Strategy const& strategy, Result& result) + { + result.distance + = strategy.apply_point_point(point1, point2); + // The projected point makes not really sense in point-point. + // We just assign one on the other + geometry::convert(point1, result.closest_point); + } +}; + +template +struct point_range +{ + template + static inline void apply(Point const& point, Range const& range, Strategy const& strategy, Result& result) + { + // This should not occur (see exception on empty input below) + if (boost::begin(range) == boost::end(range)) + { + return; + } + + // line of one point: same case as point-point + typedef typename boost::range_const_iterator::type iterator_type; + iterator_type it = boost::begin(range); + iterator_type prev = it++; + if (it == boost::end(range)) + { + point_point::type>::apply(point, *prev, strategy, result); + return; + } + + // Initialize with first segment + strategy.apply(point, *prev, *it, result); + + // Check other segments + for(prev = it++; it != boost::end(range); prev = it++) + { + Result other; + strategy.apply(point, *prev, *it, other); + if (other.distance < result.distance) + { + result = other; + } + } + } +}; + + + + +}} // namespace detail::closest_point +#endif // DOXYGEN_NO_DETAIL + + +#ifndef DOXYGEN_NO_DISPATCH +namespace dispatch +{ + +template +< + typename Geometry1, typename Geometry2, + typename Tag1 = typename tag::type, + typename Tag2 = typename tag::type, + bool Reverse = reverse_dispatch::type::value +> +struct closest_point : not_implemented +{ +}; + + +template +< + typename Geometry1, typename Geometry2, + typename Tag1, typename Tag2 +> +struct closest_point +{ + template + static inline void apply(Geometry1 const& geometry1, Geometry2 const& geometry2, + Strategy const& strategy, Result& result) + { + // Reversed version just calls dispatch with reversed arguments + closest_point + < + Geometry2, Geometry1, Tag2, Tag1, false + >::apply(geometry2, geometry1, strategy, result); + } + +}; + + +template +struct closest_point + < + Point1, Point2, + point_tag, point_tag, false + > : public detail::closest_point::point_point +{}; + + +template +struct closest_point + < + Point, Segment, + point_tag, segment_tag, + false + > +{ + template + static inline void apply(Point const& point, Segment const& segment, + Strategy const& strategy, Result& result) + { + + typename point_type::type p[2]; + geometry::detail::assign_point_from_index<0>(segment, p[0]); + geometry::detail::assign_point_from_index<1>(segment, p[1]); + + strategy.apply(point, p[0], p[1], result); + } +}; + + +template +struct closest_point + : detail::closest_point::point_range +{}; + +// +template +struct closest_point + < + Point, Linestring, + point_tag, linestring_tag, + false + > + : detail::closest_point::point_range +{}; + + +} // namespace dispatch +#endif // DOXYGEN_NO_DISPATCH + + + +template +inline void closest_point(Geometry1 const& geometry1, Geometry2 const& geometry2, Result& result) +{ + concepts::check(); + concepts::check(); + concepts::check(); + + assert_dimension_equal(); + assert_dimension_equal(); + + detail::throw_on_empty_input(geometry1); + detail::throw_on_empty_input(geometry2); + + strategy::distance::calculate_closest_point<> info_strategy; + + + dispatch::closest_point + < + Geometry1, + Geometry2 + >::apply(geometry1, geometry2, info_strategy, result); +} + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_ALGORITHMS_CLOSEST_POINT_HPP diff --git a/deps/boost/geometry/extensions/strategies/cartesian/closest_point.hpp b/deps/boost/geometry/extensions/strategies/cartesian/closest_point.hpp new file mode 100644 index 000000000..e92812e3a --- /dev/null +++ b/deps/boost/geometry/extensions/strategies/cartesian/closest_point.hpp @@ -0,0 +1,164 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2007-2013 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2008-2013 Bruno Lalande, Paris, France. +// Copyright (c) 2009-2013 Mateusz Loskot, London, UK. +// Copyright (c) 2013 Adam Wulkiewicz, Lodz, Poland. +// Copyright (c) 2017 Artem Pavlenko, UK + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_STRATEGY_CARTESIAN_CLOSEST_POINT_HPP +#define BOOST_GEOMETRY_STRATEGY_CARTESIAN_CLOSEST_POINT_HPP + +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include + + +namespace boost { namespace geometry +{ + +template +struct closest_point_result +{ + typedef Point point_type; + typedef typename default_distance_result::type distance_type; + distance_type distance; + Point closest_point; + + inline closest_point_result() + : distance(distance_type()) + {} +}; + + +namespace strategy { namespace distance +{ + +template +< + typename CalculationType = void, + typename Strategy = pythagoras +> +struct calculate_closest_point +{ +public : + // The three typedefs below are necessary to calculate distances + // from segments defined in integer coordinates. + + // Integer coordinates can still result in FP distances. + // There is a division, which must be represented in FP. + // So promote. + template + struct calculation_type + : promote_floating_point + < + typename strategy::distance::services::return_type + < + Strategy, + Point, + PointOfSegment + >::type + > + {}; + + +public : + + // Helper function + template + inline typename calculation_type::type + apply_point_point(Point1 const& p1, Point2 const& p2) const + { + Strategy point_point_strategy; + boost::ignore_unused_variable_warning(point_point_strategy); + return point_point_strategy.apply(p1, p2); + } + + template + inline void apply(Point const& p, + PointOfSegment const& p1, PointOfSegment const& p2, + Result& result) const + { + assert_dimension_equal(); + + typedef typename calculation_type::type calculation_type; + + //// A projected point of points in Integer coordinates must be able to be + //// represented in FP. + typedef model::point + < + calculation_type, + dimension::value, + typename coordinate_system::type + > fp_point_type; + + // For convenience + typedef fp_point_type fp_vector_type; + fp_vector_type v, w; + + geometry::convert(p2, v); + geometry::convert(p, w); + subtract_point(v, p1); + subtract_point(w, p1); + + calculation_type const zero = calculation_type(); + + calculation_type const c1 = dot_product(w, v); + calculation_type const c2 = dot_product(v, v); + + if (c1 < zero) + { + result.distance = apply_point_point(p, p1); + result.closest_point = p1; + } + else if(c1 > c2) + { + result.distance = apply_point_point(p, p2); + result.closest_point = p2; + } + else + { + Strategy point_point_strategy; + boost::ignore_unused_variable_warning(point_point_strategy); + + if (geometry::math::equals(c2, zero)) + { + geometry::convert(p1, result.closest_point); + result.distance = apply_point_point(p, p1); + return; + } + + calculation_type const b = c1 / c2; + geometry::convert(p1, result.closest_point); + multiply_value(v, b); + add_point(result.closest_point, v); + result.distance = apply_point_point(p, result.closest_point); + } + } +}; + +}} // namespace strategy::distance + + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGY_CARTESIAN_CLOSEST_POINT_HPP diff --git a/include/mapnik/geometry/closest_point.hpp b/include/mapnik/geometry/closest_point.hpp new file mode 100644 index 000000000..4dd7a2989 --- /dev/null +++ b/include/mapnik/geometry/closest_point.hpp @@ -0,0 +1,144 @@ +/***************************************************************************** + * + * This file is part of Mapnik (c++ mapping toolkit) + * + * Copyright (C) 2017 Artem Pavlenko + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + * + *****************************************************************************/ + +#ifndef MAPNIK_GEOMETRY_CLOSEST_POINT_HPP +#define MAPNIK_GEOMETRY_CLOSEST_POINT_HPP + +#include + +#if BOOST_VERSION >= 106200 + +#include +#include +#include +#include +#include +#include + +namespace mapnik { namespace geometry { + +namespace detail { + +template +struct closest_point +{ + using coordinate_type = T; + using result_type = boost::geometry::closest_point_result>; + + closest_point(mapnik::geometry::point const& pt) + : pt_(pt) {} + + result_type operator() (mapnik::geometry::geometry_empty const&) const + { + return result_type(); // FIXME: consider std::optional + } + + result_type operator() (mapnik::geometry::point const& pt) const + { + result_type info; + boost::geometry::closest_point(pt_ ,pt, info); + return info; + } + + result_type operator() (mapnik::geometry::line_string const& line) const + { + result_type info; + boost::geometry::closest_point(pt_ ,line, info); + return info; + } + + result_type operator() (mapnik::geometry::polygon const& poly) const + { + result_type info; + if (boost::geometry::within(pt_, poly)) + { + info.closest_point = pt_; + info.distance = 0.0; + return info; + } + bool first = true; + for (auto const& ring : poly) + { + result_type ring_info; + boost::geometry::closest_point(pt_ ,ring, ring_info); + if (first) + { + first = false; + info = std::move(ring_info); + } + else if (ring_info.distance < info.distance) + { + info = std::move(ring_info); + } + } + return info; + } + + + + // Multi* + GeometryCollection + result_type operator() (mapnik::geometry::geometry const& geom) const + { + return mapnik::util::apply_visitor(*this, geom); + } + + template + result_type operator() (MultiGeometry const& multi_geom) const + { + result_type info; + bool first = true; + for (auto const& geom : multi_geom) + { + if (first) + { + first = false; + info = std::move(operator()(geom)); + } + else + { + auto sub_info = operator()(geom); + if (sub_info.distance < info.distance) + { + info = std::move(sub_info); + } + } + } + return info; + } + + mapnik::geometry::point pt_; +}; + +} + +template +inline typename detail::closest_point::result_type +closest_point(T1 const& geom, mapnik::geometry::point const& pt) +{ + return detail::closest_point(pt)(geom); +} + +}} + + +#endif // +#endif // MAPNIK_GEOMETRY_CLOSEST_POINT_HPP diff --git a/test/unit/geometry/closest_point.cpp b/test/unit/geometry/closest_point.cpp new file mode 100644 index 000000000..2d1fd8427 --- /dev/null +++ b/test/unit/geometry/closest_point.cpp @@ -0,0 +1,46 @@ +#include "catch.hpp" + +#include + +TEST_CASE("geometry closest point") { + +SECTION("point") { + + mapnik::geometry::point pt(0, 0); + mapnik::geometry::point geom(3.0, 4.0); + auto result = mapnik::geometry::closest_point(pt, geom); + REQUIRE(result.closest_point == geom); + REQUIRE(result.distance == 5.0); +} + +SECTION("linestring") { + + //mapnik::geometry::line_string line; + //line.emplace_back(0, 0); + //line.emplace_back(25, 25); + //line.emplace_back(50, 50); + //mapnik::geometry::point centroid; + //REQUIRE(mapnik::geometry::centroid(line, centroid)); + //REQUIRE(centroid.x == 25); + //REQUIRE(centroid.y == 25); +} + + +SECTION("polygon") { + + //mapnik::geometry::polygon poly; + //mapnik::geometry::linear_ring ring; + //ring.emplace_back(0, 0); + //ring.emplace_back(1, 0); + //ring.emplace_back(1, 1); + //ring.emplace_back(0, 1); + //ring.emplace_back(0, 0); + //poly.push_back(std::move(ring)); + + //mapnik::geometry::point centroid; + //REQUIRE(mapnik::geometry::centroid(poly, centroid)); + ///REQUIRE(centroid.x == 0.5); + //REQUIRE(centroid.y == 0.5); +} + +}