Merge pull request #155 from awulkiew/fix/centroid

Fix/centroid
This commit is contained in:
Adam Wulkiewicz
2014-10-18 17:35:29 +02:00
4 changed files with 266 additions and 22 deletions

View File

@@ -5,6 +5,11 @@
// Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
// Copyright (c) 2014 Adam Wulkiewicz, Lodz, Poland.
// This file was modified by Oracle on 2014.
// Modifications copyright (c) 2014 Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
@@ -19,7 +24,6 @@
#include <cstddef>
#include <boost/range.hpp>
#include <boost/type_traits/remove_reference.hpp>
#include <boost/variant/static_visitor.hpp>
#include <boost/variant/apply_visitor.hpp>
#include <boost/variant/variant_fwd.hpp>
@@ -36,6 +40,7 @@
#include <boost/geometry/geometries/concepts/check.hpp>
#include <boost/geometry/algorithms/assign.hpp>
#include <boost/geometry/algorithms/detail/interior_iterator.hpp>
#include <boost/geometry/algorithms/convert.hpp>
#include <boost/geometry/algorithms/not_implemented.hpp>
@@ -50,6 +55,8 @@
#include <boost/geometry/algorithms/num_points.hpp>
#include <boost/geometry/multi/algorithms/num_points.hpp>
#include <boost/geometry/algorithms/detail/centroid/translating_transformer.hpp>
namespace boost { namespace geometry
{
@@ -188,17 +195,19 @@ inline bool range_ok(Range const& range, Point& centroid)
return true;
}
/*!
\brief Calculate the centroid of a ring.
\brief Calculate the centroid of a Ring or a Linestring.
*/
template <closure_selector Closure>
struct centroid_range_state
{
template<typename Ring, typename Strategy>
template<typename Ring, typename PointTransformer, typename Strategy>
static inline void apply(Ring const& ring,
Strategy const& strategy, typename Strategy::state_type& state)
PointTransformer const& transformer,
Strategy const& strategy,
typename Strategy::state_type& state)
{
typedef typename geometry::point_type<Ring const>::type point_type;
typedef typename closeable_view<Ring const, Closure>::type view_type;
typedef typename boost::range_iterator<view_type const>::type iterator_type;
@@ -207,11 +216,19 @@ struct centroid_range_state
iterator_type it = boost::begin(view);
iterator_type end = boost::end(view);
for (iterator_type previous = it++;
it != end;
++previous, ++it)
typename PointTransformer::result_type
previous_pt = transformer.apply(*it);
for ( ++it ; it != end ; ++it)
{
strategy.apply(*previous, *it, state);
typename PointTransformer::result_type
pt = transformer.apply(*it);
strategy.apply(static_cast<point_type const&>(previous_pt),
static_cast<point_type const&>(pt),
state);
previous_pt = pt;
}
}
};
@@ -221,13 +238,20 @@ struct centroid_range
{
template<typename Range, typename Point, typename Strategy>
static inline void apply(Range const& range, Point& centroid,
Strategy const& strategy)
Strategy const& strategy)
{
if (range_ok(range, centroid))
{
// prepare translation transformer
translating_transformer<Range> transformer(*boost::begin(range));
typename Strategy::state_type state;
centroid_range_state<Closure>::apply(range, strategy, state);
centroid_range_state<Closure>::apply(range, transformer,
strategy, state);
strategy.result(state, centroid);
// translate the result back
transformer.apply_reverse(centroid);
}
}
};
@@ -240,14 +264,16 @@ struct centroid_range
*/
struct centroid_polygon_state
{
template<typename Polygon, typename Strategy>
template<typename Polygon, typename PointTransformer, typename Strategy>
static inline void apply(Polygon const& poly,
Strategy const& strategy, typename Strategy::state_type& state)
PointTransformer const& transformer,
Strategy const& strategy,
typename Strategy::state_type& state)
{
typedef typename ring_type<Polygon>::type ring_type;
typedef centroid_range_state<geometry::closure<ring_type>::value> per_ring;
per_ring::apply(exterior_ring(poly), strategy, state);
per_ring::apply(exterior_ring(poly), transformer, strategy, state);
typename interior_return_type<Polygon const>::type
rings = interior_rings(poly);
@@ -255,7 +281,7 @@ struct centroid_polygon_state
for (typename detail::interior_iterator<Polygon const>::type
it = boost::begin(rings); it != boost::end(rings); ++it)
{
per_ring::apply(*it, strategy, state);
per_ring::apply(*it, transformer, strategy, state);
}
}
};
@@ -268,9 +294,16 @@ struct centroid_polygon
{
if (range_ok(exterior_ring(poly), centroid))
{
// prepare translation transformer
translating_transformer<Polygon>
transformer(*boost::begin(exterior_ring(poly)));
typename Strategy::state_type state;
centroid_polygon_state::apply(poly, strategy, state);
centroid_polygon_state::apply(poly, transformer, strategy, state);
strategy.result(state, centroid);
// translate the result back
transformer.apply_reverse(centroid);
}
}
};
@@ -282,11 +315,14 @@ struct centroid_polygon
*/
struct centroid_multi_point_state
{
template <typename Point, typename Strategy>
template <typename Point, typename PointTransformer, typename Strategy>
static inline void apply(Point const& point,
Strategy const& strategy, typename Strategy::state_type& state)
PointTransformer const& transformer,
Strategy const& strategy,
typename Strategy::state_type& state)
{
strategy.apply(point, state);
strategy.apply(static_cast<Point const&>(transformer.apply(point)),
state);
}
};
@@ -303,8 +339,9 @@ template <typename Policy>
struct centroid_multi
{
template <typename Multi, typename Point, typename Strategy>
static inline void apply(Multi const& multi, Point& centroid,
Strategy const& strategy)
static inline void apply(Multi const& multi,
Point& centroid,
Strategy const& strategy)
{
#if ! defined(BOOST_GEOMETRY_CENTROID_NO_THROW)
// If there is nothing in any of the ranges, it is not possible
@@ -315,6 +352,9 @@ struct centroid_multi
}
#endif
// prepare translation transformer
translating_transformer<Multi> transformer(multi);
typename Strategy::state_type state;
for (typename boost::range_iterator<Multi const>::type
@@ -322,9 +362,12 @@ struct centroid_multi
it != boost::end(multi);
++it)
{
Policy::apply(*it, strategy, state);
Policy::apply(*it, transformer, strategy, state);
}
Strategy::result(state, centroid);
// translate the result back
transformer.apply_reverse(centroid);
}
};

View File

@@ -0,0 +1,119 @@
// Boost.Geometry (aka GGL, Generic Geometry Library)
// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
// Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
// Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
// Copyright (c) 2014 Adam Wulkiewicz, Lodz, Poland.
// This file was modified by Oracle on 2014.
// Modifications copyright (c) 2014 Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
// 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_DETAIL_CENTROID_TRANSLATING_TRANSFORMER_HPP
#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_CENTROID_TRANSLATING_TRANSFORMER_HPP
#include <cstddef>
#include <boost/core/addressof.hpp>
#include <boost/core/ref.hpp>
#include <boost/geometry/core/cs.hpp>
#include <boost/geometry/core/tag_cast.hpp>
#include <boost/geometry/core/tags.hpp>
#include <boost/geometry/core/point_type.hpp>
#include <boost/geometry/arithmetic/arithmetic.hpp>
#include <boost/geometry/iterators/point_iterator.hpp>
namespace boost { namespace geometry
{
#ifndef DOXYGEN_NO_DETAIL
namespace detail { namespace centroid
{
// NOTE: There is no need to translate in other coordinate systems than
// cartesian. But if it was needed then one should translate using
// CS-specific technique, e.g. in spherical/geographic a translation
// vector should contain coordinates being multiplies of 2PI or 360 deg.
template <typename Geometry,
typename CastedTag = typename tag_cast
<
typename tag<Geometry>::type,
areal_tag
>::type,
typename CSTag = typename cs_tag<Geometry>::type>
struct translating_transformer
{
typedef typename geometry::point_type<Geometry>::type point_type;
typedef boost::reference_wrapper<point_type const> result_type;
explicit translating_transformer(Geometry const&) {}
explicit translating_transformer(point_type const&) {}
result_type apply(point_type const& pt) const
{
return result_type(pt);
}
template <typename ResPt>
void apply_reverse(ResPt &) const {}
};
// Specialization for Areal Geometries in cartesian CS
template <typename Geometry>
struct translating_transformer<Geometry, areal_tag, cartesian_tag>
{
typedef typename geometry::point_type<Geometry>::type point_type;
typedef point_type result_type;
explicit translating_transformer(Geometry const& geom)
: m_origin(NULL)
{
geometry::point_iterator<Geometry const>
pt_it = geometry::points_begin(geom);
if ( pt_it != geometry::points_end(geom) )
{
m_origin = boost::addressof(*pt_it);
}
}
explicit translating_transformer(point_type const& origin)
: m_origin(boost::addressof(origin))
{}
result_type apply(point_type const& pt) const
{
point_type res = pt;
if ( m_origin )
geometry::subtract_point(res, *m_origin);
return res;
}
template <typename ResPt>
void apply_reverse(ResPt & res_pt) const
{
if ( m_origin )
geometry::add_point(res_pt, *m_origin);
}
const point_type * m_origin;
};
}} // namespace detail::centroid
#endif // DOXYGEN_NO_DETAIL
}} // namespace boost::geometry
#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_CENTROID_TRANSLATING_TRANSFORMER_HPP

View File

@@ -5,6 +5,11 @@
// Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
// Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
// This file was modified by Oracle on 2014.
// Modifications copyright (c) 2014 Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
@@ -129,6 +134,34 @@ void test_large_integers()
BOOST_CHECK_EQUAL(bg::get<1>(int_centroid), bg::get<1>(double_centroid_as_int));
}
//#include <to_svg.hpp>
void test_large_doubles()
{
typedef bg::model::point<double, 2, bg::cs::cartesian> point;
point pt_far, pt_near;
bg::model::polygon<point> poly_far, poly_near;
// related to ticket #10643
bg::read_wkt("POLYGON((1074699.93 703064.65, 1074703.90 703064.58, 1074704.53 703061.40, 1074702.10 703054.62, 1074699.93 703064.65))", poly_far);
bg::read_wkt("POLYGON((699.93 64.65, 703.90 64.58, 704.53 61.40, 702.10 54.62, 699.93 64.65))", poly_near);
bg::centroid(poly_far, pt_far);
bg::centroid(poly_near, pt_near);
BOOST_CHECK(bg::within(pt_far, poly_far));
BOOST_CHECK(bg::within(pt_near, poly_near));
point pt_near_moved;
bg::set<0>(pt_near_moved, bg::get<0>(pt_near) + 1074000.0);
bg::set<1>(pt_near_moved, bg::get<1>(pt_near) + 703000.0);
//geom_to_svg(poly_far, pt_far, "far.svg");
//geom_to_svg(poly_near, pt_near, "near.svg");
double d = bg::distance(pt_far, pt_near_moved);
BOOST_CHECK(d < 0.1);
}
int test_main(int, char* [])
{
@@ -149,6 +182,9 @@ int test_main(int, char* [])
// The test currently fails in release mode. TODO: fix this
test_large_integers();
#endif
test_large_doubles();
test_exceptions<bg::model::d2::point_xy<double> >();
return 0;

View File

@@ -187,6 +187,52 @@ inline void turns_to_svg(Turns const& turns, Mapper & mapper, bool /*enrich*/ =
}
}
template <typename G1, typename P>
inline void geom_to_svg(G1 const& g1, bg::svg_mapper<P> & mapper)
{
mapper.add(g1);
mapper.map(g1, "fill-opacity:0.5;fill:rgb(153,204,0);"
"stroke:rgb(153,204,0);stroke-width:3");
}
template <typename G1, typename G2, typename P>
inline void geom_to_svg(G1 const& g1, G2 const& g2, bg::svg_mapper<P> & mapper)
{
mapper.add(g1);
mapper.add(g2);
mapper.map(g1, "fill-opacity:0.5;fill:rgb(153,204,0);"
"stroke:rgb(153,204,0);stroke-width:3");
mapper.map(g2, "fill-opacity:0.3;fill:rgb(51,51,153);"
"stroke:rgb(51,51,153);stroke-width:3");
}
template <typename G1, typename G2>
inline void geom_to_svg(G1 const& g1, std::string const& filename)
{
namespace bg = boost::geometry;
typedef typename bg::point_type<G1>::type mapper_point_type;
std::ofstream svg(filename.c_str(), std::ios::trunc);
bg::svg_mapper<mapper_point_type> mapper(svg, 500, 500);
geom_to_svg(g1, mapper);
}
template <typename G1, typename G2>
inline void geom_to_svg(G1 const& g1, G2 const& g2, std::string const& filename)
{
namespace bg = boost::geometry;
typedef typename bg::point_type<G1>::type mapper_point_type;
std::ofstream svg(filename.c_str(), std::ios::trunc);
bg::svg_mapper<mapper_point_type> mapper(svg, 500, 500);
geom_to_svg(g1, g2, mapper);
}
struct to_svg_assign_policy
: bg::detail::overlay::assign_null_policy
{