From 3b51da2d6f4c68ffc9d8c3e34ba28464836c6a41 Mon Sep 17 00:00:00 2001 From: Barend Gehrels Date: Fri, 1 Nov 2013 19:49:06 +0000 Subject: [PATCH] [geometry] added the (not yet finished) distance_info to extensions [SVN r86538] --- extensions/test/algorithms/Jamfile.v2 | 1 + extensions/test/algorithms/distance_info.cpp | 132 ++++++++++ .../extensions/algorithms/distance_info.hpp | 232 ++++++++++++++++++ .../strategies/cartesian/distance_info.hpp | 180 ++++++++++++++ 4 files changed, 545 insertions(+) create mode 100644 extensions/test/algorithms/distance_info.cpp create mode 100644 include/boost/geometry/extensions/algorithms/distance_info.hpp create mode 100644 include/boost/geometry/extensions/strategies/cartesian/distance_info.hpp diff --git a/extensions/test/algorithms/Jamfile.v2 b/extensions/test/algorithms/Jamfile.v2 index a5570a225..494a8f7be 100644 --- a/extensions/test/algorithms/Jamfile.v2 +++ b/extensions/test/algorithms/Jamfile.v2 @@ -11,6 +11,7 @@ test-suite boost-geometry-extensions-algorithms : [ run dissolve.cpp ] + [ run distance_info.cpp ] [ run connect.cpp ] [ run offset.cpp ] [ run midpoints.cpp ] diff --git a/extensions/test/algorithms/distance_info.cpp b/extensions/test/algorithms/distance_info.cpp new file mode 100644 index 000000000..ddbe2ef84 --- /dev/null +++ b/extensions/test/algorithms/distance_info.cpp @@ -0,0 +1,132 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// Unit Test + +// 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. + +// 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) + +#include + +#include + +#include +#include +#include +#include +#include +#include + + +template +void check_distance_info(Result const& result, + std::string const& expected_pp, + bool expected_on_segment, + double expected_projected_distance, + double expected_real_distance, + double expected_fraction) +{ + if (! expected_pp.empty()) + { + std::ostringstream out; + out << bg::wkt(result.projected_point1); + std::string wkt_projected = out.str(); + + BOOST_CHECK_EQUAL(wkt_projected, expected_pp); + } + BOOST_CHECK_EQUAL(result.on_segment, expected_on_segment); + BOOST_CHECK_CLOSE(result.fraction1, expected_fraction, 0.001); + BOOST_CHECK_CLOSE(result.projected_distance1, expected_projected_distance, 0.001); + BOOST_CHECK_CLOSE(result.real_distance, expected_real_distance, 0.001); +} + +template +void test_distance_info(Geometry1 const& geometry1, Geometry2 const& geometry2, + std::string const& expected_pp, + bool expected_on_segment, + double expected_projected_distance, + double expected_real_distance, + double expected_fraction) +{ + typename bg::distance_info_result::type> result, reversed_result; + bg::distance_info(geometry1, geometry2, result); + check_distance_info(result, + expected_pp, expected_on_segment, + expected_projected_distance, expected_real_distance, + expected_fraction); + + // Check reversed version too. + std::string reversed_expected_pp = expected_pp; + if (boost::is_same::type, bg::point_tag>::value + && boost::is_same::type, bg::point_tag>::value + ) + { + // For point-point, we cannot check projected-point again, it is also the other one. + reversed_expected_pp.clear(); + } + bg::distance_info(geometry2, geometry1, reversed_result); + check_distance_info(reversed_result, + reversed_expected_pp, + expected_on_segment, + expected_projected_distance, expected_real_distance, + expected_fraction); +} + +template +void test_distance_info(std::string const& wkt, std::string const& wkt_point, + std::string const& expected_pp, + bool expected_on_segment, + double expected_projected_distance, + double expected_real_distance, + double expected_fraction) +{ + Geometry1 geometry1; + typename bg::point_type::type point; + bg::read_wkt(wkt, geometry1); + bg::read_wkt(wkt_point, point); + + test_distance_info(geometry1, point, expected_pp, expected_on_segment, + expected_projected_distance, expected_real_distance, + expected_fraction); +} + +template +void test_2d() +{ + test_distance_info >("LINESTRING(2 0,4 0)", "POINT(3 2)", "POINT(3 0)", true, 2.0, 2.0, 0.5); + test_distance_info >("LINESTRING(2 0,4 0)", "POINT(2 0)", "POINT(2 0)", true, 0.0, 0.0, 0.0); + test_distance_info >("LINESTRING(2 0,4 0)", "POINT(4 0)", "POINT(4 0)", true, 0.0, 0.0, 1.0); + test_distance_info >("LINESTRING(2 0,4 0)", "POINT(5 2)", "POINT(5 0)", false, 2.0, sqrt(5.0), 1.5); + test_distance_info >("LINESTRING(2 0,4 0)", "POINT(0 2)", "POINT(0 0)", false, 2.0, sqrt(8.0), -1.0); + + // Degenerated segment + test_distance_info >("LINESTRING(2 0,2 0)", "POINT(4 0)", "POINT(2 0)", false, 2.0, 2.0, 0.0); + + // Linestring + test_distance_info >("LINESTRING(2 0,4 0)", "POINT(3 2)", "POINT(3 0)", true, 2.0, 2.0, 0.5); + + + // Point-point + test_distance_info

("Point(1 1)", "POINT(2 2)", "POINT(2 2)", false, sqrt(2.0), sqrt(2.0), 0.0); +} + +template +void test_3d() +{ + test_distance_info >("LINESTRING(0 0 0,5 5 5)", "POINT(2 3 4)", "POINT(3 3 3)", true, sqrt(2.0), sqrt(2.0), 0.6); +} + +int test_main(int, char* []) +{ + test_2d >(); + test_2d >(); + test_2d >(); + + test_3d >(); + + return 0; +} diff --git a/include/boost/geometry/extensions/algorithms/distance_info.hpp b/include/boost/geometry/extensions/algorithms/distance_info.hpp new file mode 100644 index 000000000..c174ab1c0 --- /dev/null +++ b/include/boost/geometry/extensions/algorithms/distance_info.hpp @@ -0,0 +1,232 @@ +// 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. + +// 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_DISTANCE_INFO_HPP +#define BOOST_GEOMETRY_ALGORITHMS_DISTANCE_INFO_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 distance_info { + + +template +struct point_point +{ + template + static inline void apply(Point1 const& point1, Point2 const& point2, Strategy const& strategy, Result& result) + { + result.real_distance + = result.projected_distance1 + = result.projected_distance2 + = 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.projected_point2); + geometry::convert(point2, result.projected_point1); + } +}; + +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.real_distance < result.real_distance) + { + result = other; + } + } + } +}; + + + + +}} // namespace detail::distance_info +#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 distance_info : not_implemented +{ +}; + + +template +< + typename Geometry1, typename Geometry2, + typename Tag1, typename Tag2 +> +struct distance_info +{ + template + static inline void apply(Geometry1 const& geometry1, Geometry2 const& geometry2, + Strategy const& strategy, Result& result) + { + // Reversed version just calls dispatch with reversed arguments + distance_info + < + Geometry2, Geometry1, Tag2, Tag1, false + >::apply(geometry2, geometry1, strategy, result); + } + +}; + + +template +struct distance_info + < + Point1, Point2, + point_tag, point_tag, false + > : public detail::distance_info::point_point +{}; + + +template +struct distance_info + < + 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 +//< +// typename Point, typename Ring, +// typename Point +//> +//struct distance_info +// < +// point_tag, ring_tag, +// Point, Ring, +// Point +// > +// : detail::distance_info::point_range +//{}; +// +// +template +struct distance_info + < + Point, Linestring, + point_tag, linestring_tag, + false + > + : detail::distance_info::point_range +{}; + + +} // namespace dispatch +#endif // DOXYGEN_NO_DISPATCH + + + +template +inline void distance_info(Geometry1 const& geometry1, Geometry2 const& geometry2, Result& result) +{ + concept::check(); + concept::check(); + concept::check(); + + assert_dimension_equal(); + assert_dimension_equal(); + + detail::throw_on_empty_input(geometry1); + detail::throw_on_empty_input(geometry2); + + strategy::distance::calculate_distance_info<> info_strategy; + + + dispatch::distance_info + < + Geometry1, + Geometry2 + >::apply(geometry1, geometry2, info_strategy, result); +} + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_ALGORITHMS_DISTANCE_INFO_HPP diff --git a/include/boost/geometry/extensions/strategies/cartesian/distance_info.hpp b/include/boost/geometry/extensions/strategies/cartesian/distance_info.hpp new file mode 100644 index 000000000..152fdd45e --- /dev/null +++ b/include/boost/geometry/extensions/strategies/cartesian/distance_info.hpp @@ -0,0 +1,180 @@ +// 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. + +// 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_DISTANCE_INFO_HPP +#define BOOST_GEOMETRY_STRATEGY_CARTESIAN_DISTANCE_INFO_HPP + +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include + + +namespace boost { namespace geometry +{ + +template +struct distance_info_result +{ + typedef Point point_type; + typedef typename default_distance_result::type distance_type; + + bool on_segment; + bool within_geometry; + distance_type real_distance; + Point projected_point1; // A on B + Point projected_point2; // B on A + distance_type projected_distance1; + distance_type projected_distance2; + distance_type fraction1; + distance_type fraction2; + segment_identifier seg_id1; + segment_identifier seg_id2; + + inline distance_info_result() + : on_segment(false) + , within_geometry(false) + , fraction1(distance_type()) + , fraction2(distance_type()) + , real_distance(distance_type()) + , projected_distance1(distance_type()) + , projected_distance2(distance_type()) + {} +}; + + +namespace strategy { namespace distance +{ + +template +< + typename CalculationType = void, + typename Strategy = pythagoras +> +struct calculate_distance_info +{ +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; + + + + // For source-code-comments, see "cartesian/distance_distance_info.hpp" + 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); + + result.on_segment = c1 >= zero && c1 <= c2; + + Strategy point_point_strategy; + boost::ignore_unused_variable_warning(point_point_strategy); + + if (geometry::math::equals(c2, zero)) + { + geometry::convert(p1, result.projected_point1); + result.fraction1 = 0.0; + result.on_segment = false; + result.projected_distance1 = result.real_distance = apply_point_point(p, p1); + return; + } + + calculation_type const b = c1 / c2; + result.fraction1 = b; + + geometry::convert(p1, result.projected_point1); + multiply_value(v, b); + add_point(result.projected_point1, v); + result.projected_distance1 = apply_point_point(p, result.projected_point1); + result.real_distance + = c1 < zero ? apply_point_point(p, p1) + : c1 > c2 ? apply_point_point(p, p2) + : result.projected_distance1; + } +}; + +}} // namespace strategy::distance + + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGY_CARTESIAN_DISTANCE_INFO_HPP