[Boost-commit] svn:boost r86554 - in trunk: boost/geometry/algorithms boost/geometry/algorithms/detail libs/geometry/test/algorithms libs/geometry/test/algorithms/overlay

Subject: [Boost-commit] svn:boost r86554 - in trunk: boost/geometry/algorithms boost/geometry/algorithms/detail libs/geometry/test/algorithms libs/geometry/test/algorithms/overlay
From: barend.gehrels_at_[hidden]
Date: 2013-11-03 16:00:35

Author: barendgehrels
Date: 2013-11-03 16:00:34 EST (Sun, 03 Nov 2013)
New Revision: 86554
URL: http://svn.boost.org/trac/boost/changeset/86554

Log:
[geometry] added point_on_surface, developed last summer to SVN

Added:
   trunk/boost/geometry/algorithms/detail/extreme_points.hpp   (contents, props changed)
   trunk/boost/geometry/algorithms/point_on_surface.hpp   (contents, props changed)
   trunk/libs/geometry/test/algorithms/point_on_surface.cpp   (contents, props changed)
Text files modified: 
   trunk/boost/geometry/algorithms/detail/extreme_points.hpp     |   490 ++++++++++++++++++++++++++++++++++++++++
   trunk/boost/geometry/algorithms/point_on_surface.hpp          |   349 ++++++++++++++++++++++++++++            
   trunk/libs/geometry/test/algorithms/Jamfile.v2                |     1                                         
   trunk/libs/geometry/test/algorithms/overlay/overlay_cases.hpp |    10                                         
   trunk/libs/geometry/test/algorithms/point_on_surface.cpp      |   297 ++++++++++++++++++++++++                
   5 files changed, 1147 insertions(+), 0 deletions(-)

Added: trunk/boost/geometry/algorithms/detail/extreme_points.hpp
==============================================================================
--- /dev/null	00:00:00 1970	(empty, because file is newly added)
+++ trunk/boost/geometry/algorithms/detail/extreme_points.hpp	2013-11-03 16:00:34 EST (Sun, 03 Nov 2013)	(r86554)
@@ -0,0 +1,490 @@
+// 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_DETAIL_EXTREME_POINTS_HPP
+#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_EXTREME_POINTS_HPP
+
+
+#include <cstddef>
+
+#include <boost/range.hpp>
+#include <boost/typeof/typeof.hpp>
+
+#include <boost/geometry/core/cs.hpp>
+#include <boost/geometry/core/point_type.hpp>
+#include <boost/geometry/core/ring_type.hpp>
+
+#include <boost/geometry/geometries/concepts/check.hpp>
+#include <boost/geometry/iterators/ever_circling_iterator.hpp>
+    
+#include <boost/geometry/strategies/side.hpp>
+
+#include <boost/geometry/util/math.hpp>
+
+
+namespace boost { namespace geometry
+{
+
+
+#ifndef DOXYGEN_NO_DETAIL
+namespace detail { namespace extreme_points
+{
+
+template <std::size_t Dimension>
+struct compare
+{
+    template <typename Point>
+    inline bool operator()(Point const& lhs, Point const& rhs)
+    {
+        return geometry::get<Dimension>(lhs) < geometry::get<Dimension>(rhs);
+    }
+};
+
+
+template <std::size_t Dimension, typename PointType, typename CoordinateType>
+inline void move_along_vector(PointType& point, PointType const& extreme, CoordinateType const& base_value)
+{
+    // Moves a point along the vector (point, extreme) in the direction of the extreme  point
+    // This adapts the possibly uneven legs of the triangle (or trapezium-like shape) 
+    //      _____extreme            _____
+    //     /     \                 /     \
+    //    /base   \         =>    /       \ point
+    //             \ point
+    //
+    // For so-called intruders, it can be used to adapt both legs to the level of "base"
+    // For the base, it can be used to adapt both legs to the level of the max-value of the intruders
+    // If there are 2 or more extreme values, use the one close to 'point' to have a correct vector
+
+    CoordinateType const value = geometry::get<Dimension>(point);
+    //if (geometry::math::equals(value, base_value))
+    if (value >= base_value)
+    {
+        return;
+    }
+
+    PointType vector = point;
+    subtract_point(vector, extreme);
+
+    CoordinateType const diff = geometry::get<Dimension>(vector);
+
+    // diff should never be zero
+    // because of the way our triangle/trapezium is build. 
+    // We just return if it would be the case.
+    if (geometry::math::equals(diff, 0))
+    {
+        return;
+    }
+
+    CoordinateType const base_diff = base_value - geometry::get<Dimension>(extreme);
+
+    multiply_value(vector, base_diff / diff);
+
+    // The real move:
+    point = extreme;
+    add_point(point, vector);
+}
+
+
+template <std::size_t Dimension, typename Range, typename CoordinateType>
+inline void move_along_vector(Range& range, CoordinateType const& base_value)
+{
+    if (range.size() >= 3)
+    {
+        move_along_vector<Dimension>(range.front(), *(range.begin() + 1), base_value);
+        move_along_vector<Dimension>(range.back(), *(range.rbegin() + 1), base_value);
+    }
+}
+
+
+template<typename Ring, std::size_t Dimension>
+struct extreme_points_on_ring
+{
+
+    typedef typename geometry::coordinate_type<Ring>::type coordinate_type;
+    typedef typename boost::range_iterator<Ring const>::type range_iterator;
+    typedef typename geometry::point_type<Ring>::type point_type;
+
+    typedef typename geometry::strategy::side::services::default_strategy
+        <
+            typename geometry::cs_tag<point_type>::type
+        >::type side_strategy;
+
+
+    template <typename CirclingIterator, typename Points>
+    static inline bool extend(CirclingIterator& it,
+            std::size_t n, 
+            coordinate_type max_coordinate_value, 
+            Points& points, int direction)
+    {
+        std::size_t safe_index = 0;
+        do
+        {
+            it += direction;
+
+            points.push_back(*it);
+
+            if (safe_index++ >= n)
+            {
+                // E.g.: ring is completely horizontal or vertical (= invalid, but we don't want to have an infinite loop)
+                return false;
+            }
+        } while (geometry::math::equals(geometry::get<Dimension>(*it), max_coordinate_value));
+
+        return true;
+    }
+
+    // Overload without adding to poinst
+    template <typename CirclingIterator>
+    static inline bool extend(CirclingIterator& it,
+            std::size_t n, 
+            coordinate_type max_coordinate_value, 
+            int direction)
+    {
+        std::size_t safe_index = 0;
+        do
+        {
+            it += direction;
+
+            if (safe_index++ >= n)
+            {
+                // E.g.: ring is completely horizontal or vertical (= invalid, but we don't want to have an infinite loop)
+                return false;
+            }
+        } while (geometry::math::equals(geometry::get<Dimension>(*it), max_coordinate_value));
+
+        return true;
+    }
+
+    template <typename CirclingIterator>
+    static inline bool extent_both_sides(Ring const& ring, 
+            point_type extreme,
+            CirclingIterator& left,
+            CirclingIterator& right)
+    {
+        std::size_t const n = boost::size(ring);
+        coordinate_type const max_coordinate_value = geometry::get<Dimension>(extreme);
+
+        if (! extend(left, n, max_coordinate_value, -1))
+        {
+            return false;
+        }
+        if (! extend(right, n, max_coordinate_value, +1))
+        {
+            return false;
+        }
+
+        return true;
+    }
+
+    template <typename Collection, typename CirclingIterator>
+    static inline bool collect(Ring const& ring, 
+            point_type extreme,
+            Collection& points, 
+            CirclingIterator& left,
+            CirclingIterator& right)
+    {
+        std::size_t const n = boost::size(ring);
+        coordinate_type const max_coordinate_value = geometry::get<Dimension>(extreme);
+
+        // Collects first left, which is reversed (if more than one point) then adds the top itself, then right
+        if (! extend(left, n, max_coordinate_value, points, -1))
+        {
+            return false;
+        }
+        std::reverse(points.begin(), points.end());
+        points.push_back(extreme);
+        if (! extend(right, n, max_coordinate_value, points, +1))
+        {
+            return false;
+        }
+
+        return true;
+    }
+
+    template <typename Extremes, typename Intruders, typename CirclingIterator>
+    static inline void get_intruders(Ring const& ring, CirclingIterator left, CirclingIterator right,
+            Extremes const& extremes,
+            Intruders& intruders)
+    {
+        if (boost::size(extremes) < 3)
+        {
+            return;
+        }
+        coordinate_type const min_value = geometry::get<Dimension>(*std::min_element(boost::begin(extremes), boost::end(extremes), compare<Dimension>()));
+            
+        // Also select left/right (if Dimension=1)
+        coordinate_type const other_min = geometry::get<1 - Dimension>(*std::min_element(boost::begin(extremes), boost::end(extremes), compare<1 - Dimension>()));
+        coordinate_type const other_max = geometry::get<1 - Dimension>(*std::max_element(boost::begin(extremes), boost::end(extremes), compare<1 - Dimension>()));
+        
+        std::size_t defensive_check_index = 0; // in case we skip over left/right check, collect modifies right too
+        std::size_t const n = boost::size(ring);
+        while (left != right && defensive_check_index < n)
+        {
+            coordinate_type const coordinate = geometry::get<Dimension>(*right);
+            coordinate_type const other_coordinate = geometry::get<1 - Dimension>(*right);
+            if (coordinate > min_value && other_coordinate > other_min && other_coordinate < other_max)
+            {
+                int const first_side = side_strategy::apply(*right, extremes.front(), *(extremes.begin() + 1));
+                int const last_side = side_strategy::apply(*right, *(extremes.rbegin() + 1), extremes.back());
+
+                // If not lying left from any of the extemes side
+                if (first_side != 1 && last_side != 1)
+                {
+                    //std::cout << "first " << first_side << " last " << last_side << std::endl;
+
+                    // we start at this intrusion until it is handled, and don't affect our initial left iterator
+                    CirclingIterator left_intrusion_it = right;
+                    point_type local_top = *right;
+                    typename boost::range_value<Intruders>::type intruder;
+                    collect(ring, *right, intruder, left_intrusion_it, right);
+
+                    // Also moves these to base-level, makes sorting possible which can be done in case of self-tangencies
+                    // (we might postpone this action, it is often not necessary. However it is not time-consuming)
+                    move_along_vector<Dimension>(intruder, min_value);
+                    intruders.push_back(intruder);
+                    --right;
+                }
+            }
+            ++right;
+            defensive_check_index++;
+        }
+    }
+
+    template <typename Extremes, typename Intruders>
+    static inline void get_intruders(Ring const& ring,
+            Extremes const& extremes,
+            Intruders& intruders)
+    {
+        std::size_t const n = boost::size(ring);
+        if (n >= 3)
+        {
+            geometry::ever_circling_range_iterator<Ring const> left(ring);
+            geometry::ever_circling_range_iterator<Ring const> right(ring);
+            ++right;
+
+            get_intruders(ring, left, right, extremes, intruders);
+        }
+    }
+
+    template <typename Iterator>
+    static inline bool right_turn(Ring const& ring, Iterator it)
+    {
+        int const index = std::distance(boost::begin(ring), it);
+        geometry::ever_circling_range_iterator<Ring const> left(ring);
+        geometry::ever_circling_range_iterator<Ring const> right(ring);
+        left += index;
+        right += index;
+
+        if (! extent_both_sides(ring, *it, left, right))
+        {
+            return false;
+        }
+
+        int const first_side = side_strategy::apply(*(right - 1), *right, *left);
+        int const last_side = side_strategy::apply(*left, *(left + 1), *right);
+
+//std::cout << "Candidate at " << geometry::wkt(*it) << " first=" << first_side << " last=" << last_side << std::endl;
+
+        // Turn should not be left (actually, it should be right because extent removes horizontal/collinear cases)
+        return first_side != 1 && last_side != 1;
+    }
+
+
+    // Gets the extreme segments (top point plus neighbouring points), plus intruders, if any, on the same ring
+    template <typename Extremes, typename Intruders>
+    static inline bool apply(Ring const& ring, Extremes& extremes, Intruders& intruders)
+    {
+        std::size_t const n = boost::size(ring);
+        if (n < 3)
+        {
+            return false;
+        }
+
+        // Get all maxima, usually one. In case of self-tangencies, or self-crossings,
+        // the max might be is not valid. A valid max should make a right turn
+        range_iterator max_it = boost::begin(ring);
+        compare<Dimension> smaller;
+        for (range_iterator it = max_it + 1; it != boost::end(ring); ++it)
+        {
+            if (smaller(*max_it, *it) && right_turn(ring, it))
+            {
+                max_it = it;
+            }
+        }
+
+        if (max_it == boost::end(ring))
+        {
+            return false;
+        }
+       
+        int const index = std::distance(boost::begin(ring), max_it);
+//std::cout << "Extreme point lies at " << index << " having " << geometry::wkt(*max_it) << std::endl;
+
+        geometry::ever_circling_range_iterator<Ring const> left(ring);
+        geometry::ever_circling_range_iterator<Ring const> right(ring);
+        left += index;
+        right += index;
+
+        // Collect all points (often 3) in a temporary vector
+        std::vector<point_type> points;
+        points.reserve(3);
+        if (! collect(ring, *max_it, points, left, right))
+        {
+            return false;
+        }
+
+//std::cout << "Built vector of " << points.size() << std::endl;
+
+        coordinate_type const front_value = geometry::get<Dimension>(points.front());
+        coordinate_type const back_value = geometry::get<Dimension>(points.back());
+        coordinate_type const base_value = (std::max)(front_value, back_value);
+        if (front_value < back_value)
+        {
+            move_along_vector<Dimension>(points.front(), *(points.begin() + 1), base_value);
+        }
+        else
+        {
+            move_along_vector<Dimension>(points.back(), *(points.rbegin() + 1), base_value);
+        }
+
+        std::copy(points.begin(), points.end(), std::back_inserter(extremes));
+
+        get_intruders(ring, left, right, extremes, intruders);
+
+        return true;
+    }
+};
+
+
+
+
+
+}} // namespace detail::extreme_points
+#endif // DOXYGEN_NO_DETAIL
+
+
+#ifndef DOXYGEN_NO_DISPATCH
+namespace dispatch
+{
+
+
+template
+<
+    typename Geometry,
+    std::size_t Dimension,
+    typename GeometryTag = typename tag<Geometry>::type
+>
+struct extreme_points
+{};
+
+
+template<typename Ring, std::size_t Dimension>
+struct extreme_points<Ring, Dimension, ring_tag>
+    : detail::extreme_points::extreme_points_on_ring<Ring, Dimension>
+{};
+
+
+template<typename Polygon, std::size_t Dimension>
+struct extreme_points<Polygon, Dimension, polygon_tag>
+{
+    template <typename Extremes, typename Intruders>
+    static inline bool apply(Polygon const& polygon, Extremes& extremes, Intruders& intruders)
+    {
+        typedef typename geometry::ring_type<Polygon>::type ring_type;
+        typedef detail::extreme_points::extreme_points_on_ring
+            <
+                ring_type, Dimension
+            > ring_implementation;
+
+        if (! ring_implementation::apply(geometry::exterior_ring(polygon), extremes, intruders))
+        {
+            return false;
+        }
+
+        // For a polygon, its interior rings can contain intruders
+        typename interior_return_type<Polygon const>::type rings
+                    = interior_rings(polygon);
+        for (BOOST_AUTO_TPL(it, boost::begin(rings));
+            it != boost::end(rings);
+            ++it)
+        {
+            ring_implementation::get_intruders(*it, extremes,  intruders);
+        }
+
+        return true;
+    }
+};
+
+template<typename Box>
+struct extreme_points<Box, 1, box_tag>
+{
+    template <typename Extremes, typename Intruders>
+    static inline bool apply(Box const& box, Extremes& extremes, Intruders& )
+    {
+        extremes.resize(4);
+        geometry::detail::assign_box_corners_oriented<false>(box, extremes);
+        // ll,ul,ur,lr, contains too exactly the right info
+        return true;
+    }
+};
+
+template<typename Box>
+struct extreme_points<Box, 0, box_tag>
+{
+    template <typename Extremes, typename Intruders>
+    static inline bool apply(Box const& box, Extremes& extremes, Intruders& )
+    {
+        extremes.resize(4);
+        geometry::detail::assign_box_corners_oriented<false>(box, extremes);
+        // ll,ul,ur,lr, rotate one to start with UL and end with LL
+        std::rotate(extremes.begin(), extremes.begin() + 1, extremes.end());
+        return true;
+    }
+};
+
+} // namespace dispatch
+#endif // DOXYGEN_NO_DISPATCH
+
+
+/*!
+\brief Returns extreme points (for Edge=1 in dimension 1, so the top, 
+       for Edge=0 in dimension 0, the right side)
+\note We could specify a strategy (less/greater) to get bottom/left side too. However, until now we don't need that.
+ */
+template <std::size_t Edge, typename Geometry, typename Extremes, typename Intruders>
+inline bool extreme_points(Geometry const& geometry, Extremes& extremes, Intruders& intruders)
+{
+    concept::check<Geometry const>();
+
+    // Extremes is not required to follow a geometry concept (but it should support an output iterator), 
+    // but its elements should fulfil the point-concept
+    concept::check<typename boost::range_value<Extremes>::type>();
+
+    // Intruders should contain collections which value type is point-concept
+    // Extremes might be anything (supporting an output iterator), but its elements should fulfil the point-concept
+    concept::check
+        <
+            typename boost::range_value
+                <
+                    typename boost::range_value<Intruders>::type
+                >::type
+            const
+        >();
+
+    return dispatch::extreme_points<Geometry, Edge>::apply(geometry, extremes, intruders);
+}
+
+
+
+}} // namespace boost::geometry
+
+
+#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_EXTREME_POINTS_HPP

Added: trunk/boost/geometry/algorithms/point_on_surface.hpp
==============================================================================
--- /dev/null	00:00:00 1970	(empty, because file is newly added)
+++ trunk/boost/geometry/algorithms/point_on_surface.hpp	2013-11-03 16:00:34 EST (Sun, 03 Nov 2013)	(r86554)
@@ -0,0 +1,349 @@
+// 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_POINT_ON_SURFACE_HPP
+#define BOOST_GEOMETRY_ALGORITHMS_POINT_ON_SURFACE_HPP
+
+
+#include <cstddef>
+
+#include <numeric>
+
+#include <boost/concept_check.hpp>
+#include <boost/range.hpp>
+
+#include <boost/geometry/core/point_type.hpp>
+#include <boost/geometry/core/ring_type.hpp>
+
+#include <boost/geometry/geometries/concepts/check.hpp>
+
+#include <boost/geometry/algorithms/detail/extreme_points.hpp>
+
+#include <boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp>
+
+
+namespace boost { namespace geometry
+{
+
+
+#ifndef DOXYGEN_NO_DETAIL
+namespace detail { namespace point_on_surface
+{
+
+template <typename CoordinateType, int Dimension>
+struct specific_coordinate_first
+{
+    CoordinateType const m_value_to_be_first;
+
+    inline specific_coordinate_first(CoordinateType value_to_be_skipped)
+        : m_value_to_be_first(value_to_be_skipped)
+    {}
+
+    template <typename Point>
+    inline bool operator()(Point const& lhs, Point const& rhs)
+    {
+        CoordinateType const lh = geometry::get<Dimension>(lhs);
+        CoordinateType const rh = geometry::get<Dimension>(rhs);
+
+        // If both lhs and rhs equal m_value_to_be_first, 
+        // we should handle conform if lh < rh = FALSE
+        // The first condition meets that, keep it first
+        if (geometry::math::equals(rh, m_value_to_be_first))
+        {
+            // Handle conform lh < -INF, which is always false
+            return false;
+        }
+
+        if (geometry::math::equals(lh, m_value_to_be_first))
+        {
+            // Handle conform -INF < rh, which is always true
+            return true;
+        }
+
+        return lh < rh;
+    }
+};
+
+template <int Dimension, typename Collection, typename Value, typename Predicate>
+inline bool max_value(Collection const& collection, Value& the_max, Predicate const& predicate)
+{
+    bool first = true;
+    for (typename Collection::const_iterator it = collection.begin(); it != collection.end(); ++it)
+    {
+        if (! it->empty())
+        {
+            Value the_value = geometry::get<Dimension>(*std::max_element(it->begin(), it->end(), predicate));
+            if (first || the_value > the_max)
+            {
+                the_max = the_value;
+                first = false;
+            }
+        }
+    }
+    return ! first;
+}
+
+
+template <int Dimension, typename Value>
+struct select_below
+{
+    Value m_value;
+    inline select_below(Value const& v)
+        : m_value(v)
+    {}
+
+    template <typename Intruder>
+    inline bool operator()(Intruder const& intruder) const
+    {
+        if (intruder.empty())
+        {
+            return true;
+        }
+        Value max = geometry::get<Dimension>(*std::max_element(intruder.begin(), intruder.end(), detail::extreme_points::compare<Dimension>()));
+        return geometry::math::equals(max, m_value) || max < m_value;
+    }
+};
+
+template <int Dimension, typename Value>
+struct adapt_base
+{
+    Value m_value;
+    inline adapt_base(Value const& v)
+        : m_value(v)
+    {}
+
+    template <typename Intruder>
+    inline void operator()(Intruder& intruder) const
+    {
+        if (intruder.size() >= 3)
+        {
+            detail::extreme_points::move_along_vector<Dimension>(intruder, m_value);
+        }
+    }
+};
+
+template <int Dimension, typename Value>
+struct min_of_intruder
+{
+    template <typename Intruder>
+    inline bool operator()(Intruder const& lhs, Intruder const& rhs) const
+    {
+        Value lhs_min = geometry::get<Dimension>(*std::min_element(lhs.begin(), lhs.end(), detail::extreme_points::compare<Dimension>()));
+        Value rhs_min = geometry::get<Dimension>(*std::min_element(rhs.begin(), rhs.end(), detail::extreme_points::compare<Dimension>()));
+        return lhs_min < rhs_min;
+    }
+};
+
+template <typename Point, typename Segments>
+inline void calculate_centroid(Point& point, Segments const& segments)
+{
+    if (segments.size() == 3)
+    {
+        // In almost all cases, the segments do have 3 values. In that case we use another
+        // centroid calculation, which should be slightly faster 
+        // and is more precise (case #geos_1_test_overlay => normal centroid is outside! TODO)
+
+        typedef typename geometry::coordinate_type<Point>::type coordinate_type;
+        coordinate_type const three = 3.0;
+        coordinate_type const sx = geometry::get<0>(segments[0]) + geometry::get<0>(segments[1]) + geometry::get<0>(segments[2]);
+        coordinate_type const sy = geometry::get<1>(segments[0]) + geometry::get<1>(segments[1]) + geometry::get<1>(segments[2]);
+        geometry::set<0>(point, sx / three);
+        geometry::set<1>(point, sy / three);
+        return;
+    }
+
+    // For 4 or more segments, we use normal centroid calculation
+
+    // Specify centroid, it should have "areal_tag" to have correct calculation
+    typedef typename strategy::centroid::services::default_strategy
+        <
+            typename cs_tag<Point>::type,
+            areal_tag,
+            dimension<Point>::type::value,
+            Point,
+            Point
+        >::type strategy_type;
+
+    strategy_type strategy;
+
+
+    // Ignore warning (because using static method sometimes) on strategy
+    boost::ignore_unused_variable_warning(strategy);
+
+
+    typename strategy_type::state_type state;
+
+    typedef typename boost::range_iterator<Segments const>::type iterator_type;
+
+    iterator_type begin = boost::begin(segments);
+    iterator_type end = boost::end(segments);
+    iterator_type it = begin;
+    iterator_type previous = it++;
+    for (; it != end; ++previous, ++it)
+    {
+        strategy.apply(*previous, *it, state);
+    }
+    // Close it explicitly:
+    strategy.apply(*previous, *begin, state);
+
+    strategy.result(state, point);
+}
+
+
+
+template <int Dimension, typename Extremes, typename Intruders, typename CoordinateType>
+inline void replace_extremes_for_self_tangencies(Extremes& extremes, Intruders& intruders, CoordinateType const& max_intruder)
+{
+    // This function handles self-tangencies.
+    // Self-tangencies use, as usual, the major part of code...
+
+    //        ___ e
+    //       /|\ \          
+    //      / | \ \        
+    //     /  |  \ \       
+    //    /   |   \ \
+    //   / /\ |    \ \
+    //     i2    i1
+
+    // The picture above shows the extreme (outside, "e") and two intruders ("i1","i2")
+    // Assume that "i1" is self-tangent with the extreme, in one point at the top
+    // Now the "penultimate" value is searched, this is is the top of i2
+    // Then everything including and below (this is "i2" here) is removed
+    // Then the base of "i1" and of "e" is adapted to this penultimate value
+    // It then looks like:
+
+    //      b ___ e
+    //       /|\ \          
+    //      / | \ \        
+    //     /  |  \ \       
+    //    /   |   \ \
+    //   a    c i1
+
+    // Then intruders (here "i1" but there may be more) are sorted from left to right
+    // Finally points "a","b" and "c" (in this order) are selected as a new triangle.
+    // This triangle will have a centroid which is inside (assumed that intruders left segment 
+    // is not equal to extremes left segment, but that polygon would be invalid)
+
+    // Find highest non-self tangent intrusion, if any
+    CoordinateType penultimate_value;
+    specific_coordinate_first<CoordinateType, Dimension> pu_compare(max_intruder);
+    if (max_value<Dimension>(intruders, penultimate_value, pu_compare))
+    {
+        // Throw away all intrusions <= this value, and of the kept one set this as base.
+        select_below<Dimension, CoordinateType> predicate(penultimate_value);
+        intruders.erase
+            (
+                std::remove_if(boost::begin(intruders), boost::end(intruders), predicate),
+                boost::end(intruders)
+            );
+        adapt_base<Dimension, CoordinateType> fe_predicate(penultimate_value);
+        // Sort from left to right (or bottom to top if Dimension=0)
+        std::for_each(boost::begin(intruders), boost::end(intruders), fe_predicate);
+
+        // Also adapt base of extremes
+        detail::extreme_points::move_along_vector<Dimension>(extremes, penultimate_value);
+    }
+    // Then sort in 1-Dim. Take first to calc centroid.
+    std::sort(boost::begin(intruders), boost::end(intruders), min_of_intruder<1 - Dimension, CoordinateType>());
+
+    Extremes triangle;
+    triangle.reserve(3);
+
+    // Make a triangle of first two points of extremes (the ramp, from left to right), and last point of first intruder (which goes from right to left)
+    std::copy(extremes.begin(), extremes.begin() + 2, std::back_inserter(triangle));
+    triangle.push_back(intruders.front().back());
+
+    // (alternatively we could use the last two points of extremes, and first point of last intruder...):
+    //// ALTERNATIVE: std::copy(extremes.rbegin(), extremes.rbegin() + 2, std::back_inserter(triangle));
+    //// ALTERNATIVE: triangle.push_back(intruders.back().front());
+    
+    // Now replace extremes with this smaller subset, a triangle, such that centroid calculation will result in a point inside
+    extremes = triangle;
+}
+
+template <int Dimension, typename Geometry, typename Point>
+inline bool calculate_point_on_surface(Geometry const& geometry, Point& point)
+{
+    typedef typename geometry::point_type<Geometry>::type point_type;
+    typedef typename geometry::coordinate_type<Geometry>::type coordinate_type;
+    std::vector<point_type> extremes;
+
+    typedef std::vector<std::vector<point_type> > intruders_type;
+    intruders_type intruders;
+    geometry::extreme_points<Dimension>(geometry, extremes, intruders);
+
+    if (extremes.size() < 3)
+    {
+        return false;
+    }
+
+    // If there are intruders, find the max.
+    if (! intruders.empty())
+    {
+        coordinate_type max_intruder;
+        detail::extreme_points::compare<Dimension> compare;
+        if (max_value<Dimension>(intruders, max_intruder, compare))
+        {
+            coordinate_type max_extreme = geometry::get<Dimension>(*std::max_element(extremes.begin(), extremes.end(), detail::extreme_points::compare<Dimension>()));
+            if (max_extreme > max_intruder)
+            {
+                detail::extreme_points::move_along_vector<Dimension>(extremes, max_intruder);
+            }
+            else
+            {
+                replace_extremes_for_self_tangencies<Dimension>(extremes, intruders, max_intruder);
+            }
+        }
+    }
+
+    // Now calculate the centroid of the (possibly adapted) extremes
+    calculate_centroid(point, extremes);
+
+    return true;
+}
+
+}} // namespace detail::point_on_surface
+#endif // DOXYGEN_NO_DETAIL
+
+
+/*!
+\brief Returns point guaranteed to lie on the surface
+\tparam Geometry geometry type. This also defines the type of the output point
+\param point to assign
+\param geometry geometry to take point from
+ */
+template <typename Geometry, typename Point>
+inline void point_on_surface(Geometry const& geometry, Point& point)
+{
+    concept::check<Point>();
+    concept::check<Geometry const>();
+
+    // First try in Y-direction (which should always succeed for valid polygons)
+    if (! detail::point_on_surface::calculate_point_on_surface<1>(geometry, point))
+    {
+        // For invalid polygons, we might try X-direction
+        detail::point_on_surface::calculate_point_on_surface<0>(geometry, point);
+    }
+}
+
+
+template<typename Geometry>
+inline typename geometry::point_type<Geometry>::type return_point_on_surface(Geometry const& geometry)
+{
+    typename geometry::point_type<Geometry>::type result;
+    point_on_surface(geometry, result);
+    return result;
+}
+
+}} // namespace boost::geometry
+
+
+#endif // BOOST_GEOMETRY_ALGORITHMS_POINT_ON_SURFACE_HPP

Modified: trunk/libs/geometry/test/algorithms/Jamfile.v2
==============================================================================
--- trunk/libs/geometry/test/algorithms/Jamfile.v2	Sun Nov  3 15:46:56 2013	(r86553)
+++ trunk/libs/geometry/test/algorithms/Jamfile.v2	2013-11-03 16:00:34 EST (Sun, 03 Nov 2013)	(r86554)
@@ -33,6 +33,7 @@
     [ run make.cpp ]
     [ run overlaps.cpp ]
     [ run perimeter.cpp ]
+    [ run point_on_surface.cpp ]
     [ run remove_spikes.cpp ]
     [ run reverse.cpp ]
     [ run simplify.cpp ]

Modified: trunk/libs/geometry/test/algorithms/overlay/overlay_cases.hpp
==============================================================================
--- trunk/libs/geometry/test/algorithms/overlay/overlay_cases.hpp	Sun Nov  3 15:46:56 2013	(r86553)
+++ trunk/libs/geometry/test/algorithms/overlay/overlay_cases.hpp	2013-11-03 16:00:34 EST (Sun, 03 Nov 2013)	(r86554)
@@ -575,6 +575,16 @@
     "POLYGON((2.76143527 1.093332171 , 2.076887131 1.822299719 , 4.263789177 3.875944376 , 4.948337555 3.146976948 , 2.76143527 1.093332171))"
     };
 
+static std::string ticket_8180[2] =
+    {
+    "POLYGON(( -2.559375047683716 -0.615625500679016, -2.559375047683716 0.384374797344208, 7.940625190734863 0.384374588727951, 7.940625190734863 -0.615625441074371 ))",
+    "POLYGON(( 1.000000000000000 0.384374707937241, 1.000000000000000 0.000000000000000, 0.000000000000000 0.000000000000000, 0.000000000000000 0.384374737739563 ))"
+    };
+static std::string ticket_8180b[2] =
+    {
+    "POLYGON(( -2.559375047683716 -0.615625500679016, -2.559375047683716 0.384374797344208, 7.940625190734863 0.384374588727951, 7.940625190734863 -0.615625441074371 ))",
+    "POLYGON(( 1 0.384374707937241, 1 1, 0 1, 0 0.384374737739563 ))"
+    };
 
 static std::string ticket_8254[2] =
     {

Added: trunk/libs/geometry/test/algorithms/point_on_surface.cpp
==============================================================================
--- /dev/null	00:00:00 1970	(empty, because file is newly added)
+++ trunk/libs/geometry/test/algorithms/point_on_surface.cpp	2013-11-03 16:00:34 EST (Sun, 03 Nov 2013)	(r86554)
@@ -0,0 +1,297 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+// Unit Test
+
+// 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.
+
+// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
+// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
+
+// 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)
+
+#define TEST_WITH_SVG
+
+// Instead of having a separate (and nearly similar) unit test to test multipolygons,
+// we now include them here and compile them by default. Only undefining the next line
+// will avoid testing multi-geometries
+// #define BOOST_GEOMETRY_UNIT_TEST_MULTI
+
+// TODO: multis are not supported yet
+
+
+#include <geometry_test_common.hpp>
+
+// The include to test
+#include <boost/geometry/algorithms/point_on_surface.hpp>
+
+// Helper includes
+#include <boost/geometry/algorithms/correct.hpp>
+#include <boost/geometry/algorithms/within.hpp>
+#include <boost/geometry/io/wkt/wkt.hpp>
+#include <boost/geometry/strategies/strategies.hpp>
+
+#include <boost/geometry/geometries/geometries.hpp>
+#include <boost/geometry/geometries/point_xy.hpp>
+
+// Include from unit tests
+#include <algorithms/test_overlay.hpp>
+#include <algorithms/test_relate.hpp>
+#include <algorithms/overlay/overlay_cases.hpp>
+
+#if defined(BOOST_GEOMETRY_UNIT_TEST_MULTI)
+#  include <boost/geometry/multi/algorithms/correct.hpp>
+#  include <boost/geometry/multi/algorithms/within.hpp>
+#  include <boost/geometry/multi/geometries/multi_polygon.hpp>
+#  include <boost/geometry/multi/io/wkt/wkt.hpp>
+#endif
+
+
+#if defined(TEST_WITH_SVG)
+#  include <boost/geometry/io/svg/svg_mapper.hpp>
+#endif
+
+template <typename Geometry>
+void test_geometry(std::string const& case_id, std::string const& wkt, double expected_x, double expected_y)
+{
+//std::cout << case_id << std::endl;
+    typedef typename bg::point_type<Geometry>::type point_type;
+
+    Geometry geometry;
+    bg::read_wkt(wkt, geometry);
+    bg::correct(geometry);
+
+    point_type point;
+    bg::point_on_surface(geometry, point);
+
+    BOOST_CHECK_MESSAGE(bg::within(point, geometry), 
+        case_id << " generated point_on_surface (dim 1) is not within geometry");
+
+#ifdef BOOST_GEOMETRY_POINT_ON_SURFACE_COMPLETE
+    // For the algorithm we also check generation in the other dimension
+    point_type right_point;
+    bg::detail::point_on_surface::calculate_point_on_surface<0>(geometry, right_point);
+
+    BOOST_CHECK_MESSAGE(bg::within(right_point, geometry), 
+        case_id << " generated point_on_surface (dim 0) is not within geometry");
+
+    point_type returned_point = bg::return_point_on_surface(geometry);
+#endif
+
+#if defined(TEST_WITH_SVG)
+    {
+        std::ostringstream filename;
+        filename << "point_on_surface_" << case_id << "_"
+            << string_from_type<typename bg::coordinate_type<Geometry>::type>::name()
+            << ".svg";
+
+        std::ofstream svg(filename.str().c_str());
+
+        // To map the intermediate products:
+        bg::model::linestring<point_type> top_points;
+        std::vector<bg::model::linestring<point_type> > top_intruders;
+        bg::extreme_points<1>(geometry, top_points, top_intruders);
+
+#ifdef BOOST_GEOMETRY_POINT_ON_SURFACE_COMPLETE
+        bg::model::linestring<point_type> right_points;
+        std::vector<bg::model::linestring<point_type> > right_intruders;
+        bg::extreme_points<0>(geometry, right_points, right_intruders);
+#endif
+
+        bg::svg_mapper<point_type> mapper(svg, 500, 500);
+        mapper.add(geometry);
+        mapper.map(geometry, "fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153,204,0);stroke-width:1");
+
+        // Top (red/magenta)
+        mapper.map(top_points, "stroke:rgb(255,0,0);stroke-width:2");
+        BOOST_FOREACH(auto const& intruder, top_intruders)
+        {
+            mapper.map(intruder, "stroke:rgb(255,0,255);stroke-width:2");
+        }
+        mapper.map(point, "opacity:0.8;fill:rgb(255,128,0);stroke:rgb(0,0,100);stroke-width:1", 3);
+
+#ifdef BOOST_GEOMETRY_POINT_ON_SURFACE_COMPLETE
+        //// Right (blue/cyan)
+        // (mostly commented, makes the picture less clear)
+        //mapper.map(right_points, "stroke:rgb(0,0,255);stroke-width:2");
+        //BOOST_FOREACH(auto const& intruder, right_intruders)
+        //{
+        //    mapper.map(intruder, "stroke:rgb(0,255,255);stroke-width:2");
+        //}
+        mapper.map(right_point, "opacity:0.8;fill:rgb(0,128,255);stroke:rgb(0,0,100);stroke-width:1", 3);
+#endif
+    }
+#endif
+
+}
+
+template <typename Point>
+void test_all()
+{
+    typedef bg::model::polygon<Point> polygon;
+ 
+    // Specific test-cases for point-on-surface:
+    test_geometry<polygon>("ice", "polygon((5 0,0 5,1 6,2 5,4 8,6 5,7 6,8 5,9 6,10 5,5 0))", 0, 0);
+    test_geometry<polygon>("intruding", "polygon((5 0,0 5,1 6,2 5,4 8,6 5,7 6,8 5,9 6,10 5,7 2,4 7,6 1,5 0))", 0, 0);
+    test_geometry<polygon>("intruding2", "polygon((5 0,3 2,4 6,2 3,0 5,1 6,2 5,4 8,6 5,7 6,8 5,9 6,10 5,7 2,4 7,6 1,5 0))", 0, 0);
+    test_geometry<polygon>("intruding3", "polygon((5 0,3 2,3 6,2 3,0 5,1 6,3 7,2 5,4 8,6 5,5 7,7 6,8 5,9 6,10 5,7 2,5 6,6 1,5 0))", 0, 0);
+    test_geometry<polygon>("intruding4", "polygon((5 0,3 2,3 6,2 3,0 5,1 6,3 7,2 5,4 8,6 5,5 8,7 6,8 5,9 6,10 5,7 2,5 6,6 1,5 0))", 0, 0);
+    test_geometry<polygon>("intruding5", "polygon((5 0,3 2,3 6,2 3,0 5,1 6,3 8,2 5,4 8,6 5,5 8,7 6,8 5,9 6,10 5,7 2,5 6,6 1,5 0))", 0, 0);
+    test_geometry<polygon>("ice_int", "polygon((5 0,0 5,1 6,2 5,4 8,6 5,7 6,8 5,9 6,10 5,5 0),(3 4,4 3,5 4,4 7,3 4))", 0, 0);
+    test_geometry<polygon>("ice_int2", "polygon((5 0,0 5,1 6,2 5,4 8,6 5,7 6,8 5,9 6,10 5,5 0),(4 7,3 4,4 3,5 4,4 7))", 0, 0);
+    test_geometry<polygon>("ice_in3", "polygon((5 0,0 5,1 6,2 5,4 8,6 5,7 6,8 5,9 6,10 5,5 0),(3 4,5 4,5 6,3 6,3 4))", 0, 0);
+    test_geometry<polygon>("simplex_normal", simplex_normal[0], 0, 0);
+    test_geometry<polygon>("sqr", "polygon((0 0,0 5,0 10,5 10,10 10,10 5,10 0,5 0,0 0))", 0, 0);
+    test_geometry<polygon>("self_tangent", "polygon((0 0,5 10,10 0,9 0,7 4,8 0,7 0,5 10,3 0,2 0,3 4,1 0,0 0))", 0, 0);
+    test_geometry<polygon>("self_tangent2", "polygon((5 0,2 3,4 8,1.5 3.5,0 5,1 6,2 5,4 8,6 5,7 6,8 5,9 6,10 5,5 0))", 0, 0);
+    test_geometry<polygon>("self_tangent_int", "polygon((5 0,0 5,1 6,2 5,4 8,6 5,7 6,8 5,9 6,10 5,5 0),(3 4,4 3,5 4,4 8,3 4))", 0, 0);
+    test_geometry<polygon>("self_tangent_int2", "polygon((5 0,2 3,3.5 7,1.5 3.5,0 5,1 6,2 5,4 8,6 5,7 6,8 5,9 6,10 5,5 0),(3 4,4 3,5 4,4 8,3 4))", 0, 0);
+    test_geometry<polygon>("self_tangent_int3", "polygon((5 0,2 3,4 8,1.5 3.5,0 5,1 6,2 5,4 8,6 5,7 6,8 5,9 6,10 5,5 0),(3 4,4 3,5 4,4 8,3 4))", 0, 0);
+
+    test_geometry<polygon>("self_tangent_int3", "polygon((5 0,2 3,4 8,1.5 3.5,0 5,1 6,2 5,4 8,6 5,7 6,8 5,9 6,10 5,5 0),(3 4,4 3,5 4,4 8,3 4))", 0, 0);
+    test_geometry<polygon>("disjoint_simplex0", disjoint_simplex[0], 0, 0);
+    test_geometry<polygon>("disjoint_simplex1", disjoint_simplex[1], 0, 0);
+
+#if defined(BOOST_GEOMETRY_UNIT_TEST_MULTI)
+    typedef bg::model::multi_polygon<polygon> multi_polygon;
+    test_geometry<multi_polygon>("mp_simplx", "multipolygon(((0 0,0 5,5 0, 0 0),(7 1,7 6,10 1,7 1)))", 0, 0);
+#endif
+
+
+    // Overlay testcases
+    test_geometry<polygon>("buffer_mp1", buffer_mp1[0],  0, 0);
+    test_geometry<polygon>("buffer_mp2", buffer_mp2[0],  0, 0);
+    test_geometry<polygon>("buffer_rt_a", buffer_rt_a[0],  0, 0);
+    test_geometry<polygon>("buffer_rt_f", buffer_rt_f[0],  0, 0);
+    test_geometry<polygon>("buffer_rt_g", buffer_rt_g[0],  0, 0);
+    test_geometry<polygon>("buffer_rt_i", buffer_rt_i[0],  0, 0);
+    test_geometry<polygon>("buffer_rt_j", buffer_rt_j[0],  0, 0);
+    test_geometry<polygon>("buffer_rt_l", buffer_rt_l[0],  0, 0);
+    test_geometry<polygon>("buffer_rt_m1", buffer_rt_m1[0],  0, 0);
+    test_geometry<polygon>("buffer_rt_m2", buffer_rt_m2[0],  0, 0);
+    test_geometry<polygon>("buffer_rt_n", buffer_rt_n[0],  0, 0);
+    test_geometry<polygon>("buffer_rt_q", buffer_rt_q[0],  0, 0);
+    test_geometry<polygon>("buffer_rt_r", buffer_rt_r[0],  0, 0);
+    test_geometry<polygon>("buffer_rt_t", buffer_rt_t[0],  0, 0);
+    test_geometry<polygon>("case_10", case_10[0],  0, 0);
+    test_geometry<polygon>("case_11", case_11[0],  0, 0);
+    test_geometry<polygon>("case_12", case_12[0],  0, 0);
+    test_geometry<polygon>("case_13", case_13[0],  0, 0);
+    test_geometry<polygon>("case_14", case_14[0],  0, 0);
+    test_geometry<polygon>("case_15", case_15[0],  0, 0);
+    test_geometry<polygon>("case_16", case_16[0],  0, 0);
+    test_geometry<polygon>("case_17", case_17[0],  0, 0);
+    test_geometry<polygon>("case_18", case_18[0],  0, 0);
+    test_geometry<polygon>("case_19", case_19[0],  0, 0);
+    test_geometry<polygon>("case_1", case_1[0],  0, 0);
+    test_geometry<polygon>("case_20", case_20[0],  0, 0);
+    test_geometry<polygon>("case_21", case_21[0],  0, 0);
+    test_geometry<polygon>("case_22", case_22[0],  0, 0);
+    test_geometry<polygon>("case_23", case_23[0],  0, 0);
+    test_geometry<polygon>("case_24", case_24[0],  0, 0);
+    test_geometry<polygon>("case_25", case_25[0],  0, 0);
+    test_geometry<polygon>("case_26", case_26[0],  0, 0);
+    test_geometry<polygon>("case_27", case_27[0],  0, 0);
+    test_geometry<polygon>("case_28", case_28[0],  0, 0);
+    test_geometry<polygon>("case_29", case_29[0],  0, 0);
+    test_geometry<polygon>("case_2", case_2[0],  0, 0);
+    test_geometry<polygon>("case_30", case_30[0],  0, 0);
+    test_geometry<polygon>("case_31", case_31[0],  0, 0);
+    test_geometry<polygon>("case_32", case_32[0],  0, 0);
+    test_geometry<polygon>("case_33", case_33[0],  0, 0);
+    test_geometry<polygon>("case_34", case_34[0],  0, 0);
+    test_geometry<polygon>("case_35", case_35[0],  0, 0);
+    test_geometry<polygon>("case_36", case_36[0],  0, 0);
+    test_geometry<polygon>("case_37", case_37[0],  0, 0);
+    test_geometry<polygon>("case_38", case_38[0],  0, 0);
+    test_geometry<polygon>("case_39", case_39[0],  0, 0);
+    test_geometry<polygon>("case_3", case_3[0],  0, 0);
+    test_geometry<polygon>("case_40", case_40[0],  0, 0);
+    test_geometry<polygon>("case_41", case_41[0],  0, 0);
+    test_geometry<polygon>("case_42", case_42[0],  0, 0);
+    test_geometry<polygon>("case_43", case_43[0],  0, 0);
+    test_geometry<polygon>("case_44", case_44[0],  0, 0);
+    test_geometry<polygon>("case_45", case_45[0],  0, 0);
+    test_geometry<polygon>("case_46", case_46[0],  0, 0);
+    test_geometry<polygon>("case_47", case_47[0],  0, 0);
+    test_geometry<polygon>("case_49", case_49[0],  0, 0);
+    test_geometry<polygon>("case_4", case_4[0],  0, 0);
+    test_geometry<polygon>("case_50", case_50[0],  0, 0);
+    test_geometry<polygon>("case_51", case_51[0],  0, 0);
+    test_geometry<polygon>("case_52", case_52[0],  0, 0);
+    test_geometry<polygon>("case_53", case_53[0],  0, 0);
+    test_geometry<polygon>("case_54", case_54[0],  0, 0);
+    test_geometry<polygon>("case_55", case_55[0],  0, 0);
+    test_geometry<polygon>("case_56", case_56[0],  0, 0);
+    test_geometry<polygon>("case_57", case_57[0],  0, 0);
+    test_geometry<polygon>("case_58", case_58[0],  0, 0);
+    test_geometry<polygon>("case_59", case_59[0],  0, 0);
+    test_geometry<polygon>("case_5", case_5[0],  0, 0);
+    test_geometry<polygon>("case_60", case_60[0],  0, 0);
+    test_geometry<polygon>("case_61", case_61[0],  0, 0);
+    test_geometry<polygon>("case_6", case_6[0],  0, 0);
+    test_geometry<polygon>("case_70", case_70[0],  0, 0);
+    test_geometry<polygon>("case_71", case_71[0],  0, 0);
+    test_geometry<polygon>("case_72", case_72[0],  0, 0);
+    test_geometry<polygon>("case_79", case_79[0],  0, 0);
+    test_geometry<polygon>("case_7", case_7[0],  0, 0);
+    test_geometry<polygon>("case_8", case_8[0],  0, 0);
+    test_geometry<polygon>("case_9", case_9[0],  0, 0);
+    test_geometry<polygon>("case_many_situations", case_many_situations[0],  0, 0);
+    test_geometry<polygon>("ccw_case_1", ccw_case_1[0],  0, 0);
+    test_geometry<polygon>("ccw_case_9", ccw_case_9[0],  0, 0);
+    test_geometry<polygon>("collinear_opposite_left", collinear_opposite_left[0],  0, 0);
+    test_geometry<polygon>("collinear_opposite_right", collinear_opposite_right[0],  0, 0);
+    test_geometry<polygon>("collinear_opposite_straight", collinear_opposite_straight[0],  0, 0);
+    test_geometry<polygon>("collinear_overlaps", collinear_overlaps[0],  0, 0);
+    test_geometry<polygon>("dz_1", dz_1[0],  0, 0);
+    test_geometry<polygon>("dz_2", dz_2[0],  0, 0);
+    test_geometry<polygon>("dz_3", dz_3[0],  0, 0);
+    test_geometry<polygon>("dz_4", dz_4[0],  0, 0);
+    test_geometry<polygon>("geos_1_test_overlay", geos_1_test_overlay[0],  0, 0);
+    test_geometry<polygon>("geos_2", geos_2[0],  0, 0);
+    test_geometry<polygon>("geos_3", geos_3[0],  0, 0);
+    test_geometry<polygon>("geos_4", geos_4[0],  0, 0);
+    test_geometry<polygon>("ggl_list_20110306_javier", ggl_list_20110306_javier[0],  0, 0);
+    test_geometry<polygon>("ggl_list_20110307_javier", ggl_list_20110307_javier[0],  0, 0);
+    test_geometry<polygon>("ggl_list_20110627_phillip", ggl_list_20110627_phillip[0],  0, 0);
+    test_geometry<polygon>("ggl_list_20110716_enrico", ggl_list_20110716_enrico[0],  0, 0);
+    test_geometry<polygon>("ggl_list_20110820_christophe ", ggl_list_20110820_christophe [0],  0, 0);
+    test_geometry<polygon>("ggl_list_20120717_volker", ggl_list_20120717_volker[0],  0, 0);
+    test_geometry<polygon>("hv_1", hv_1[0],  0, 0);
+    test_geometry<polygon>("hv_2", hv_2[0],  0, 0);
+    test_geometry<polygon>("hv_3", hv_3[0],  0, 0);
+    test_geometry<polygon>("hv_4", hv_4[0],  0, 0);
+    test_geometry<polygon>("hv_5", hv_5[0],  0, 0);
+    test_geometry<polygon>("hv_6", hv_6[0],  0, 0);
+    test_geometry<polygon>("hv_7", hv_7[0],  0, 0);
+    test_geometry<polygon>("isovist", isovist[0],  0, 0);
+    test_geometry<polygon>("open_case_1", open_case_1[0],  0, 0);
+    test_geometry<polygon>("open_case_9", open_case_9[0],  0, 0);
+    test_geometry<polygon>("pie_16_2_15_0", pie_16_2_15_0[0],  0, 0);
+    test_geometry<polygon>("pie_16_4_12", pie_16_4_12[0],  0, 0);
+    test_geometry<polygon>("pie_20_20_7_100", pie_20_20_7_100[0],  0, 0);
+    test_geometry<polygon>("pie_23_16_16", pie_23_16_16[0],  0, 0);
+    test_geometry<polygon>("pie_23_21_12_500", pie_23_21_12_500[0],  0, 0);
+    test_geometry<polygon>("pie_23_23_3_2000", pie_23_23_3_2000[0],  0, 0);
+    test_geometry<polygon>("pie_4_13_15", pie_4_13_15[0],  0, 0);
+    test_geometry<polygon>("pie_5_12_12_0_7s", pie_5_12_12_0_7s[0],  0, 0);
+    test_geometry<polygon>("snl_1", snl_1[0],  0, 0);
+    test_geometry<polygon>("ticket_17", ticket_17[0],  0, 0);
+    test_geometry<polygon>("ticket_5103", ticket_5103[0],  0, 0);
+    test_geometry<polygon>("ticket_7462", ticket_7462[0],  0, 0);
+    test_geometry<polygon>("ticket_8180", ticket_8180[0],  0, 0);
+    test_geometry<polygon>("ticket_8180b", ticket_8180b[0],  0, 0);
+    test_geometry<polygon>("ticket_8254", ticket_8254[0],  0, 0);
+}
+
+
+
+int test_main(int, char* [])
+{
+    test_all<bg::model::d2::point_xy<double> >();
+
+    return 0;
+}