boost/geometry
2
0
Fork 0
mirror of https://github.com/boostorg/geometry.git synced 2026-02-12 12:12:10 +00:00
Code Issues Projects Releases Wiki Activity

Bugfix in assemble for multi

Browse Source 
Bugfix in get_turns in one constellation
Other preparations for dissolve/buffer

[SVN r59593]
This commit is contained in:
Barend Gehrels
2010-02-09 09:24:02 +00:00
parent 8f7e3593e9
commit 132fdd87ee
8 changed files with 536 additions and 479 deletions
Download Patch File Download Diff File Expand all files Collapse all files

19
include/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp
View File

@@ -47,6 +47,10 @@ struct base_turn_handler
{
ti.operations[0].operation = op;
ti.operations[1].operation = op;
if (op == operation_union)
{
ti.ignore = true;
}
}
// If condition, first union/second intersection, else vice versa
@@ -63,10 +67,7 @@ struct base_turn_handler
template <typename TurnInfo>
static inline void uu_else_ii(bool condition, TurnInfo& ti)
{
ti.operations[0].operation = condition
? operation_union : operation_intersection;
ti.operations[1].operation = condition
? operation_union : operation_intersection;
both(ti, condition ? operation_union : operation_intersection);
}
};
@@ -241,7 +242,7 @@ struct touch : public base_turn_handler
{
// Collinear -> lines join, continue
// (#BRL2)
if (side_pk_q2 == 0)
if (side_pk_q2 == 0 && ! block_q)
{
both(ti, operation_continue);
return;
@@ -292,7 +293,7 @@ struct touch : public base_turn_handler
else
{
// Pk at other side than Qi/Pk
int const side_qk_q = SideStrategy::apply(qi, qj, qk);
int const side_qk_q = SideStrategy::apply(qi, qj, qk);
bool const q_turns_left = side_qk_q == 1;
ti.operations[0].operation = q_turns_left
@@ -304,6 +305,12 @@ struct touch : public base_turn_handler
? operation_union
: operation_intersection;
if (ti.operations[0].operation == operation_union
&& ti.operations[1].operation == operation_union)
{
ti.ignore = true;
}
return;
}
}

274
include/boost/geometry/algorithms/detail/overlay/ring_properties.hpp
View File

@@ -8,12 +8,6 @@
#ifndef BOOST_GEOMETRY_ALG_DET_OV_RING_PROPERTIES_HPP
#define BOOST_GEOMETRY_ALG_DET_OV_RING_PROPERTIES_HPP
#define BOOST_GEOMETRY_ASSEMBLE_PARENT_DEQUE
#ifdef BOOST_GEOMETRY_ASSEMBLE_PARENT_DEQUE
#include <deque>
#else
#include <boost/array.hpp>
#endif
#include <boost/geometry/algorithms/area.hpp>
#include <boost/geometry/algorithms/envelope.hpp>
@@ -31,89 +25,132 @@ namespace boost { namespace geometry
template <typename Point>
struct ring_properties
{
struct parent
{
int c;
ring_identifier id;
typedef Point point_type;
typedef geometry::box<Point> box_type;
inline parent()
: c(0)
{}
ring_identifier ring_id;
typename area_result<Point>::type area;
int signum;
inline parent(int c_, ring_identifier i)
: c(c_), id(i)
{}
};
bool intersects;
bool produced;
#ifdef BOOST_GEOMETRY_ASSEMBLE_PARENT_DEQUE
inline void push(parent const& p)
// "Stack"/counter of non-intersecting parent rings.
// This denotes if it a negative ring should be included,
int parent_count;
// ID of the parent
ring_identifier parent_ring_id;
box_type box;
point_type point;
bool has_point;
// Default constructor (necessary for vector, but not called)
inline ring_properties()
: intersects(false)
, produced(false)
, parent_count(0)
, has_point(false)
{
parents.push_back(p);
parent_ring_id.source_index = -1;
}
inline void pop()
template <typename Geometry>
inline ring_properties(ring_identifier id, Geometry const& geometry,
bool i, bool p = false)
: ring_id(id)
, area(geometry::area(geometry))
, intersects(i)
, produced(p)
, parent_count(0)
, box(geometry::make_envelope<box_type>(geometry))
{
parents.pop_back();
has_point = geometry::point_on_border(point, geometry, true);
typedef typename coordinate_type<Geometry>::type coordinate_type;
coordinate_type zero = coordinate_type();
signum = area > zero ? 1 : area < zero ? -1 : 0;
parent_ring_id.source_index = -1;
}
inline bool parent_stack_empty() const
inline bool positive() const { return signum == 1; }
inline bool negative() const { return signum == -1; }
inline bool operator<(ring_properties<Point> const& other) const
{
return parents.empty();
// Normal sorting: in reverse order
return abs(area) > abs(other.area);
}
#else
inline void push(parent const& p)
inline ring_identifier const& id(int direction) const
{
parents[array_size++] = p;
// Return the id of ifself, or of the parent
return positive() || parent_ring_id.source_index < 0
? ring_id
: parent_ring_id;
}
inline void pop()
inline void push(ring_properties<Point> const& r,
int direction, bool dissolve)
{
array_size--;
}
inline bool parent_stack_empty() const
{
return array_size == 0;
}
if (//(r.produced || r.untouched()) &&
r.included(direction, dissolve))
{
#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
std::cout << " id.push " << r.ring_id;
#endif
parent_ring_id = r.ring_id;
}
if (! r.produced || dissolve)
{
#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
std::cout << " or.push " << r.ring_id;
#endif
parent_count++;
}
}
inline void pop(ring_properties<Point> const& r)
{
if (! r.produced && parent_count > 0)
{
#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
std::cout << " or.pop";
#endif
parent_count--;
}
}
inline bool interior_included(int direction) const
{
if (area < 0 && ! parent_stack_empty())
if (negative())
{
// Inner rings are included if the last encountered parent
// matches the operation
#ifdef BOOST_GEOMETRY_ASSEMBLE_PARENT_DEQUE
int d = parents.back().c;
if (d == 0
&& parents.size() == 2
&& parents.front().c == 0)
#else
int d = parents[array_size - 1].c;
if (d == 0
&& array_size == 2
&& parents[0].c == 0)
#endif
// Original inner rings are included if there
// are two untouched parents (union) or one (intersection);
// Produced are ones are included if there is a parent found
if (produced)
{
// It is contained by two outer rings, both of them
// are not included in the other. So they must be
// equal. Then always included, both in union and
// in intersection.
return direction == -1;
return parent_count > 0;
}
if (d == 0)
{
d = 1;
}
return d * direction == 1;
return
(direction == 1 && parent_count == 1)
|| (direction == -1 && parent_count > 1);
}
return false;
}
inline bool included(int direction) const
inline bool included(int direction, bool dissolve) const
{
if (produced)
if (produced && ! dissolve)
{
// Traversed rings are included in all operations,
// because traversal was direction-dependant.
// On dissolve, this is not the case.
return true;
}
if (intersects)
@@ -123,13 +160,19 @@ struct ring_properties
return false;
}
if (area > 0)
if (positive())
{
// Outer rings are included if they don't have parents
// (union) or have parents (intersection)
return (parent_stack_empty() ? 1 : -1) * direction == 1;
if (produced)
{
return parent_count == 0;
}
return
(direction == 1 && parent_count == 0)
|| (direction == -1 && parent_count > 0);
}
else if (area < 0 && ! parent_stack_empty())
else if (negative())
{
// Inner rings are included if the last encountered parent
// matches the operation
@@ -138,86 +181,28 @@ struct ring_properties
return false;
}
typedef Point point_type;
typedef geometry::box<Point> box_type;
ring_identifier ring_id;
typename area_result<Point>::type area;
int c;
bool intersects;
bool produced;
// Stack of non-intersecting parent rings.
// This denotes if it a negative ring should be included,
// and which is the parent.
#ifdef BOOST_GEOMETRY_ASSEMBLE_PARENT_DEQUE
std::deque<parent> parents;
#else
parent parents[3];
int array_size;
#endif
// If the parent is an intersecting ring, take
// this as the parent.
ring_identifier parent_id;
bool has_parent_id;
box_type box;
point_type point;
bool has_point;
inline ring_properties()
: c(0)
, intersects(false)
, produced(false)
, has_point(false)
, has_parent_id(false)
#ifndef BOOST_GEOMETRY_ASSEMBLE_PARENT_DEQUE
, array_size(0)
#endif
{}
template <typename Geometry>
inline ring_properties(ring_identifier id, Geometry const& geometry,
bool i, bool p = false)
: ring_id(id)
, area(geometry::area(geometry))
, c(0)
, intersects(i)
, produced(p)
, has_parent_id(false)
, box(geometry::make_envelope<box_type>(geometry))
#ifndef BOOST_GEOMETRY_ASSEMBLE_PARENT_DEQUE
, array_size(0)
#endif
inline bool untouched() const
{
has_point = geometry::point_on_border(point, geometry);
// It should be in comparisons on parent/child if:
// it is produced
// it is not produced, and not intersecting
return ! produced && ! intersects;
}
inline bool operator<(ring_properties<Point> const& other) const
{
// Normal sorting: in reverse order
return std::abs(area) > std::abs(other.area);
}
inline ring_identifier const& id() const
#if defined(BOOST_GEOMETRY_DEBUG_IDENTIFIER)
friend std::ostream& operator<<(std::ostream &os, ring_properties<Point> const& prop)
{
if (has_parent_id)
{
return parent_id;
}
return parent_stack_empty() || area > 0
? ring_id
#ifdef BOOST_GEOMETRY_ASSEMBLE_PARENT_DEQUE
: parents.back().id;
#else
: parents[array_size - 1].id;
os << "prop: " << prop.ring_id << " " << prop.area;
os << " count: " << prop.parent_count;
std::cout << " parent: " << prop.parent_ring_id;
if (prop.produced) std::cout << " produced";
if (prop.intersects) std::cout << " intersects";
if (prop.included(1, false)) std::cout << " union";
if (prop.included(-1, false)) std::cout << " intersection";
return os;
}
#endif
}
};
@@ -225,14 +210,21 @@ struct ring_properties
template<typename Prop>
struct sort_on_id_or_parent_id
{
private :
int m_direction;
public :
sort_on_id_or_parent_id(int direction)
: m_direction(direction)
{}
inline bool operator()(Prop const& left, Prop const& right) const
{
ring_identifier const& left_id = left.id();
ring_identifier const& right_id = right.id();
ring_identifier const& left_id = left.id(m_direction);
ring_identifier const& right_id = right.id(m_direction);
// If it is the same, sort on size descending
return left_id == right_id
? std::abs(left.area) > std::abs(right.area)
? abs(left.area) > abs(right.area)
: left_id < right_id;
}
};

4
include/boost/geometry/algorithms/detail/overlay/turn_info.hpp
View File

@@ -82,11 +82,15 @@ struct turn_info
Point point;
method_type method;
bool ignore;
bool rejected;
Container operations;
turn_info()
: method(method_none)
, ignore(false)
, rejected(false)
{}
};

2
include/boost/geometry/algorithms/intersection.hpp
View File

@@ -90,7 +90,7 @@ template
typename Strategy
>
struct intersection_inserter
: detail::overlay::overlay_and_assemble
: detail::overlay::overlay
<G1, G2, OutputIterator, GeometryOut, -1, Strategy>
{};

380
include/boost/geometry/algorithms/overlay/assemble.hpp
View File

@@ -21,7 +21,6 @@
#include <boost/geometry/algorithms/overlay/get_turns.hpp>
#include <boost/geometry/algorithms/overlay/enrich_intersection_points.hpp>
//#include <boost/geometry/algorithms/detail/overlay/enrich.hpp>
#include <boost/geometry/algorithms/overlay/traverse.hpp>
@@ -37,8 +36,6 @@
#include <boost/geometry/strategies/intersection.hpp>
//#define BOOST_GEOMETRY_DEBUG_ASSEMBLE
#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
# include <boost/geometry/util/write_dsv.hpp>
#endif
@@ -165,10 +162,16 @@ struct add_to_containment<ring_tag, Ring>
{
template <typename ContainmentContainer, typename Map>
static inline void apply(ContainmentContainer& container,
ring_identifier const& id, Ring const& ring, Map const& map)
ring_identifier const& id, Ring const& ring, Map const& map,
bool dissolve)
{
typedef typename range_value<ContainmentContainer>::type prop;
container.push_back(prop(id, ring, map.find(id) != map.end()));
bool found = map.find(id) != map.end();
if (! dissolve || ! found)
{
// For dissolve, do not add intersected rings
container.push_back(prop(id, ring, found));
}
}
};
@@ -177,10 +180,15 @@ struct add_to_containment<box_tag, Box>
{
template <typename ContainmentContainer, typename Map>
static inline void apply(ContainmentContainer& container,
ring_identifier const& id, Box const& box, Map const& map)
ring_identifier const& id, Box const& box, Map const& map,
bool dissolve)
{
typedef typename range_value<ContainmentContainer>::type prop;
container.push_back(prop(id, box, map.find(id) != map.end()));
bool found = map.find(id) != map.end();
if (! dissolve || ! found)
{
container.push_back(prop(id, box, found));
}
}
};
@@ -190,7 +198,8 @@ struct add_to_containment<polygon_tag, Polygon>
{
template <typename ContainmentContainer, typename Map>
static inline void apply(ContainmentContainer& container,
ring_identifier const& id, Polygon const& polygon, Map const& map)
ring_identifier const& id, Polygon const& polygon, Map const& map,
bool dissolve)
{
// Add exterior ring and interior rings
ring_identifier copy = id;
@@ -201,7 +210,7 @@ struct add_to_containment<polygon_tag, Polygon>
typename ring_type<Polygon>::type
> policy;
policy::apply(container, copy, exterior_ring(polygon), map);
policy::apply(container, copy, exterior_ring(polygon), map, dissolve);
copy.ring_index = 0;
for (typename boost::range_iterator
<
@@ -210,7 +219,7 @@ struct add_to_containment<polygon_tag, Polygon>
it != boost::end(interior_rings(polygon));
++it, ++copy.ring_index)
{
policy::apply(container, copy, *it, map);
policy::apply(container, copy, *it, map, dissolve);
}
}
};
@@ -229,19 +238,22 @@ inline void map_turns(Map& map, TurnPoints const& turn_points)
it != boost::end(turn_points);
++it, ++index)
{
int op_index = 0;
for (typename boost::range_iterator<container_type const>::type
op_it = boost::begin(it->operations);
op_it != boost::end(it->operations);
++op_it, ++op_index)
if (! it->ignore)
{
ring_identifier ring_id
(
op_it->seg_id.source_index,
op_it->seg_id.multi_index,
op_it->seg_id.ring_index
);
map[ring_id]++;
int op_index = 0;
for (typename boost::range_iterator<container_type const>::type
op_it = boost::begin(it->operations);
op_it != boost::end(it->operations);
++op_it, ++op_index)
{
ring_identifier ring_id
(
op_it->seg_id.source_index,
op_it->seg_id.multi_index,
op_it->seg_id.ring_index
);
map[ring_id]++;
}
}
}
}
@@ -259,46 +271,27 @@ struct enrich_containment
{
typedef typename boost::range_value<Container>::type item_type;
typedef typename item_type::parent parent_type;
typedef typename geometry::tag<Geometry1>::type tag1;
typedef typename geometry::tag<Geometry2>::type tag2;
typedef void tag3; // For the ring-container
static inline void assign(item_type& item1, item_type& item2)
static inline void assign(item_type& larger, item_type& smaller, int direction, bool dissolve)
{
int or1 = item1.area > 0 ? 1 : item1.area < 0 ? -1 : 0;
int or2 = item2.area > 0 ? 1 : item2.area < 0 ? -1 : 0;
if (or1 != 0 && or2 != 0)
{
if (or1 == 1 && or2 == 1)
{
item1.c = 1;
item2.c = -1;
}
typedef typename geometry::coordinate_type
<
Geometry1
>::type coordinate_type;
if (! item1.produced)
if (larger.signum != 0 && smaller.signum != 0)
{
if (larger.signum == 1)
{
// This is an original ring, used to set info about
// the parents of interiors
if (or1 == 1)
{
item2.push(parent_type(item1.c, item1.ring_id));
}
else if (or1 == -1 && ! item1.intersects)
{
item2.pop();
}
smaller.push(larger, direction, dissolve);
}
else
else if (larger.signum == -1)
{
if (or1 == 1 && or2 == -1 && ! item2.intersects)
{
// Assign this produced parent as the parent
// to this hole
item2.has_parent_id = true;
item2.parent_id = item1.ring_id;
}
smaller.pop(/*from*/larger);
}
}
}
@@ -348,7 +341,7 @@ struct enrich_containment
static inline void enrich(Selection& selection, Map& map,
Geometry1 const& geometry1, Geometry2 const& geometry2,
RingCollection const& collection,
int direction)
int direction, bool dissolve)
{
typedef typename boost::range_iterator<Selection>::type iterator;
@@ -370,12 +363,12 @@ struct enrich_containment
item_type& item2 = **it2;
if (geometry::within(item2.point, item1.box)
&& std::abs(item2.area) < std::abs(item1.area)
&& abs(item2.area) < abs(item1.area)
&& contains(item1, item2, geometry1, geometry2,
collection)
)
{
assign(item1, item2);
assign(item1, item2, direction, dissolve);
}
}
}
@@ -386,7 +379,7 @@ struct enrich_containment
static inline void divide_and_conquer(Selection& selection, Map& map,
Geometry1 const& geometry1, Geometry2 const& geometry2,
RingCollection const& collection,
int direction, Box const& box,
int direction, bool dissolve, Box const& box,
std::size_t iteration = 0, std::size_t previous_count = 0)
{
std::size_t n = boost::size(selection);
@@ -400,14 +393,15 @@ std::cout << "spatial divide n="
if (iteration > 3)
{
enrich(selection, map, geometry1, geometry2, collection, direction);
enrich(selection, map, geometry1, geometry2, collection, direction, dissolve);
return;
}
// Divide the complete box in two (alternating) halves
Box lower = box, upper = box;
typename geometry::coordinate_type<Box>::type two = 2.0;
typename geometry::coordinate_type<Box>::type mid
typedef typename geometry::coordinate_type<Box>::type coordinate_type;
coordinate_type two = 2.0;
coordinate_type mid
= (geometry::get<min_corner, Dimension>(box)
+ geometry::get<max_corner, Dimension>(box)) / two;
@@ -420,15 +414,15 @@ std::cout << "spatial divide n="
select_by_box(upper_sel, selection, upper);
divide_and_conquer<1 - Dimension>(lower_sel, map, geometry1, geometry2,
collection, direction, lower, iteration + 1, n);
collection, direction, dissolve, lower, iteration + 1, n);
divide_and_conquer<1 - Dimension>(upper_sel, map, geometry1, geometry2,
collection, direction, upper, iteration + 1, n);
collection, direction, dissolve, upper, iteration + 1, n);
}
static inline void enrich(Container& container,
Geometry1 const& geometry1, Geometry2 const& geometry2,
RingCollection const& collection,
int direction)
int direction, bool dissolve)
{
typedef typename boost::range_iterator<Container>::type iterator;
@@ -439,17 +433,26 @@ std::cout << "spatial divide n="
++it1)
{
item_type& item1 = *it1;
#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
std::cout << item1.ring_id << " area: " << item1.area << std::endl;
#endif
iterator it2 = it1;
for (it2++; it2 != boost::end(container); ++it2)
{
item_type& item2 = *it2;
if (geometry::within(item2.point, item1.box)
&& std::abs(item2.area) < std::abs(item1.area)
&& abs(item2.area) < abs(item1.area)
&& contains(item1, item2, geometry1, geometry2, collection)
)
{
#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
std::cout << " -> contains " << item2.ring_id;
#endif
n++;
assign(item1, item2);
assign(item1, item2, direction, dissolve);
#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
std::cout << std::endl;
#endif
}
}
}
@@ -461,13 +464,23 @@ std::cout << "spatial divide n="
static inline void apply(Container& container,
Geometry1 const& geometry1, Geometry2 const& geometry2,
RingCollection const& collection,
int direction, Box const& box)
int direction, bool dissolve, Box const& box)
{
if (boost::size(container) == 0)
{
return;
}
enrich(container, geometry1, geometry2, collection, direction);
enrich(container, geometry1, geometry2, collection, direction, dissolve);
#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
for (typename boost::range_iterator<Container const>::type
it = boost::begin(container);
it != boost::end(container);
++it)
{
std::cout << *it << std::endl;
}
#endif
return;
// Method using divide and conquer (does NOT work corretly!)
@@ -484,7 +497,7 @@ std::cout << "spatial divide n="
std::map<std::pair<item_type*, item_type*>, bool> map;
divide_and_conquer<1>(selection, map, geometry1, geometry2, collection,
direction, box);
direction, dissolve, box);
}
};
@@ -501,7 +514,7 @@ template
inline OutputIterator add_all_rings(Container& container,
Geometry1 const& geometry1, Geometry2 const& geometry2,
RingCollection const& collection,
int direction,
int direction, bool dissolve,
OutputIterator out)
{
typedef typename boost::range_iterator<Container>::type iterator;
@@ -520,7 +533,7 @@ inline OutputIterator add_all_rings(Container& container,
it != boost::end(container);
++it)
{
if (it->area > 0)
if (it->positive())
{
if (result_filled)
{
@@ -531,7 +544,7 @@ inline OutputIterator add_all_rings(Container& container,
// If it is an outer ring, it is included if there are no parents
// (union) or if there are parents (intersection)
if (it->included(direction))
if (it->included(direction, dissolve))
{
geometry::clear(result);
previous_id = it->ring_id;
@@ -558,8 +571,8 @@ inline OutputIterator add_all_rings(Container& container,
// If it is an interior ring, it is included if
// it's parent-id matches the id of the outputted exterior ring
if (result_filled
&& it->id() == previous_id
&& it->included(direction))
&& it->id(direction) == previous_id
&& it->included(direction, dissolve))
{
if (it->ring_id.source_index == 0)
{
@@ -589,140 +602,58 @@ inline OutputIterator add_all_rings(Container& container,
return out;
}
template <typename Container>
inline void sort_on_parent_id(Container& container, int direction)
{
typedef typename boost::range_iterator<Container>::type iterator;
typedef typename boost::range_value<Container>::type item_type;
// For negative rings produced by traversal,
// find the parent of the selected operation
// (by popping from the stack if necessary)
for (iterator it = boost::begin(container);
it != boost::end(container);
++it)
{
if (it->produced && it->area < 0)
{
while (! it->interior_included(direction)
&& ! it->parent_stack_empty())
{
it->pop();
}
}
}
// Sort container on parent-id
#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
std::cout << "assemble.properties sort on parent-id "
<< boost::size(container) << std::endl;
#endif
std::sort(boost::begin(container), boost::end(container),
sort_on_id_or_parent_id<item_type>());
}
template
<
typename GeometryOut,
typename Rings, typename Turns,
typename Geometry1, typename Geometry2,
typename OutputIterator, typename GeometryOut,
int Direction,
typename Strategy
typename OutputIterator
>
struct overlay_and_assemble
inline OutputIterator assemble(Rings const& rings, Turns& turn_points,
Geometry1 const& geometry1,
Geometry2 const& geometry2,
int direction, bool dissolve,
OutputIterator out)
{
typedef typename geometry::tag<Geometry1>::type tag1;
typedef typename geometry::tag<Geometry2>::type tag2;
typedef typename geometry::tag<GeometryOut>::type tag_out;
static inline OutputIterator apply(Geometry1 const& geometry1,
Geometry2 const& geometry2, OutputIterator out,
Strategy const& strategy)
{
if (geometry::num_points(geometry1) == 0
|| geometry::num_points(geometry2) == 0)
{
return out;
}
typedef typename geometry::tag<Geometry1>::type tag1;
typedef typename geometry::tag<Geometry2>::type tag2;
typedef typename geometry::tag<GeometryOut>::type tag_out;
typedef typename geometry::point_type<GeometryOut>::type point_type;
typedef typename boost::range_value<Rings>::type ring_type;
typedef detail::overlay::traversal_turn_info<point_type> turn_info;
typedef std::deque<turn_info> container_type;
// "Abuse" rangetype for ringtype: for polygon, it is the type of the
// exterior ring. For ring, it is the ring itself. That is what is
// wished here as well.
typedef typename geometry::range_type<GeometryOut>::type ring_type;
container_type turn_points;
#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
std::cout << "get turns" << std::endl;
#endif
boost::geometry::get_turns
<
detail::overlay::CalculateDistancePolicy
>(geometry1, geometry2, turn_points);
#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
std::cout << "enrich" << std::endl;
#endif
typename Strategy::side_strategy_type side_strategy;
geometry::enrich_intersection_points(turn_points, geometry1, geometry2,
side_strategy);
std::vector<ring_type> rings;
#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
std::cout << "traverse" << std::endl;
#endif
geometry::traverse<typename Strategy::side_strategy_type, ring_type>
(
geometry1,
geometry2,
Direction == -1
? boost::geometry::detail::overlay::operation_intersection
: boost::geometry::detail::overlay::operation_union
,
turn_points,
std::back_inserter(rings)
);
#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
std::cout << "assemble" << std::endl;
#endif
// Map intersection-points per ring-identifier to <count>
std::map<ring_identifier, int> map;
map_turns(map, turn_points);
// TODO: translate/implement this (sorry for the dutch)
//1. loop ook door originele ringen heen, map[ring] = 0
//2. verander dan map.find=.end in .find -> it=0
//3. we weten nu dat er "untouched" ringen zijn, als dat zo is,
// dan "enrich assemble" etc.kl
typedef std::vector
<
ring_properties<point_type>
> ring_properties_container_type;
ring_properties_container_type ring_properties_container;
add_to_containment
<
tag1,
Geometry1
>::apply(ring_properties_container,
ring_identifier(0,-1,-1), geometry1, map);
add_to_containment
<
tag2,
Geometry2
>::apply(ring_properties_container,
ring_identifier(1,-1,-1), geometry2, map);
ring_identifier(0, -1,-1), geometry1,
map, dissolve);
if (! dissolve)
{
add_to_containment
<
tag2,
Geometry2
>::apply(ring_properties_container,
ring_identifier(1, -1,-1), geometry2,
map, dissolve);
}
// Add all produced rings using source index 2
{
@@ -745,7 +676,7 @@ std::cout << "assemble" << std::endl;
it != boost::end(ring_properties_container);
++it)
{
if (it->area>0)
if (it->positive())
{
pos++;
}
@@ -788,15 +719,100 @@ std::cout << "assemble.enrich containment" << std::endl;
std::vector<ring_type>,
geometry::box<point_type>
>::apply(ring_properties_container,
geometry1, geometry2, rings, Direction, total);
geometry1, geometry2, rings, direction, dissolve, total);
sort_on_parent_id(ring_properties_container, Direction);
// Sort container on parent-id
#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
std::cout << "assemble.properties sort on parent-id "
<< boost::size(ring_properties_container) << std::endl;
#endif
std::sort(boost::begin(ring_properties_container),
boost::end(ring_properties_container),
sort_on_id_or_parent_id
<
ring_properties<point_type>
>(direction));
}
#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
std::cout << "assemble.add rings" << std::endl;
#endif
return add_all_rings<GeometryOut>(ring_properties_container,
geometry1, geometry2, rings, Direction, out);
geometry1, geometry2, rings, direction, dissolve, out);
}
template
<
typename Geometry1, typename Geometry2,
typename OutputIterator, typename GeometryOut,
int Direction,
typename Strategy
>
struct overlay
{
static inline OutputIterator apply(
Geometry1 const& geometry1, Geometry2 const& geometry2,
OutputIterator out,
Strategy const& strategy)
{
if (geometry::num_points(geometry1) == 0 && geometry::num_points(geometry2) == 0)
{
return out;
}
typedef typename geometry::point_type<GeometryOut>::type point_type;
typedef detail::overlay::traversal_turn_info<point_type> turn_info;
typedef std::deque<turn_info> container_type;
// "Use" rangetype for ringtype:
// for polygon, it is the type of the exterior ring.
// for ring, it is the ring itself. That is what is
// for multi-polygon, it is also the type of the ring.
typedef typename geometry::range_type<GeometryOut>::type ring_type;
container_type turn_points;
std::vector<ring_type> rings;
// If one input is empty, output the other one for a union.
// For an intersection, the intersection is empty.
if (geometry::num_points(geometry1) == 0
|| geometry::num_points(geometry2) == 0)
{
if (Direction == 1)
{
return assemble<GeometryOut>(rings, turn_points,
geometry1, geometry2, Direction, false, out);
}
return out;
}
#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
std::cout << "get turns" << std::endl;
#endif
boost::geometry::get_turns
<
detail::overlay::CalculateDistancePolicy
>(geometry1, geometry2, turn_points);
#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
std::cout << "enrich" << std::endl;
#endif
typename Strategy::side_strategy_type side_strategy;
geometry::enrich_intersection_points(turn_points, geometry1, geometry2,
side_strategy);
#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
std::cout << "traverse" << std::endl;
#endif
geometry::traverse(geometry1, geometry2,
Direction == -1
? boost::geometry::detail::overlay::operation_intersection
: boost::geometry::detail::overlay::operation_union
,
turn_points, rings);
return assemble<GeometryOut>(rings, turn_points,
geometry1, geometry2, Direction, false, out);
}
};

330
include/boost/geometry/algorithms/overlay/traverse.hpp
View File

@@ -18,29 +18,51 @@
#include <boost/geometry/geometries/concepts/check.hpp>
#ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION
#include <boost/geometry/extensions/gis/io/wkt/write_wkt.hpp>
#if defined(BOOST_GEOMETRY_DEBUG_INTERSECTION) || defined(BOOST_GEOMETRY_OVERLAY_REPORT_WKT)
#include <boost/geometry/extensions/gis/io/wkt/wkt.hpp>
#endif
#include <boost/geometry/algorithms/detail/overlay/turn_info.hpp>
namespace boost { namespace geometry
{
#ifndef DOXYGEN_NO_DETAIL
namespace detail { namespace overlay {
template <typename Turns>
inline void clear_visit_info(Turns& turns)
{
typedef typename boost::range_value<Turns>::type tp_type;
for (typename boost::range_iterator<Turns>::type
it = boost::begin(turns);
it != boost::end(turns);
++it)
{
for (typename boost::range_iterator
<
typename tp_type::container_type
>::type op_it = boost::begin(it->operations);
op_it != boost::end(it->operations);
++op_it)
{
op_it->visited.set_none();
}
}
}
template <typename Info, typename Turn>
inline void set_visited_for_contine(Info& info, Turn const& turn)
{
// On "continue", set "visited" for ALL directions
if (turn.operation == detail::overlay::operation_continue)
{
for (typename boost::range_iterator
@@ -64,19 +86,16 @@ template
typename GeometryOut,
typename G1,
typename G2,
typename TurnPoints,
typename Turns,
typename IntersectionInfo
>
inline void assign_next_ip(G1 const& g1, G2 const& g2,
TurnPoints& turn_points,
typename boost::range_iterator<TurnPoints>::type & ip,
Turns& turns,
typename boost::range_iterator<Turns>::type & ip,
GeometryOut& current_output,
IntersectionInfo & info,
segment_identifier& seg_id)
{
#ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION
//std::cout << "traverse: take: " << info << std::endl;
#endif
info.visited.set_visited();
set_visited_for_contine(*ip, info);
@@ -95,12 +114,12 @@ inline void assign_next_ip(G1 const& g1, G2 const& g2,
info.enriched.travels_to_vertex_index,
current_output);
}
ip = boost::begin(turn_points) + info.enriched.travels_to_ip_index;
ip = boost::begin(turns) + info.enriched.travels_to_ip_index;
seg_id = info.seg_id;
}
else
{
ip = boost::begin(turn_points) + info.enriched.next_ip_index;
ip = boost::begin(turns) + info.enriched.next_ip_index;
seg_id = info.seg_id;
}
*(std::back_inserter(current_output)++) = ip->point;
@@ -109,9 +128,7 @@ inline void assign_next_ip(G1 const& g1, G2 const& g2,
inline bool select_source(operation_type operation, int source1, int source2)
{
return
(operation == operation_intersection && source1 != source2)
|| (operation == operation_union && source1 == source2);
return operation == operation_intersection && source1 != source2;
}
@@ -134,7 +151,7 @@ inline bool select_next_ip(operation_type operation,
{
// In some cases there are two alternatives.
// For "ii", take the other one (alternate)
// For "uu", take the same one
// For "uu", take the same one (see above);
// For "cc", take either one, but if there is a starting one,
// take that one.
if ( (it->operation == operation_continue
@@ -157,34 +174,50 @@ inline bool select_next_ip(operation_type operation,
if (it->visited.started())
{
selected = it;
has_tp = true;
return true;
}
}
return has_tp;
}
template <typename Container>
inline void stop_gracefully(Container& container, bool& stop,
std::string const& reason)
// This backtracking approach is necessary for invalid (== self-interescting)
// geometries and the approach will be enhanced. TODO
template
<
typename Rings,
typename Turns,
typename Geometry1,
typename Geometry2
>
inline void backtrack(std::size_t size_at_start, bool& fail,
Rings& rings, typename boost::range_value<Rings>::type& ring,
Turns& turns, typename boost::range_value<Turns>::type& turn,
std::string const& reason,
Geometry1 const& geometry1,
Geometry2 const& geometry2)
{
#ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION
std::cout << "STOPPING: " << reason << std::endl;
fail = true;
// Make bad output clean
rings.resize(size_at_start);
ring.clear();
// Reject this as a starting point
turn.rejected = true;
// And clear all visit info
clear_visit_info(turns);
#ifdef BOOST_GEOMETRY_OVERLAY_REPORT_WKT
std::cout << " BT (" << reason << " )";
std::cout
<< geometry::wkt(geometry1) << std::endl
<< geometry::wkt(geometry2) << std::endl;
#endif
stop = true;
if (container.size() > 0)
{
// Two-step adding first value, without assignment
// and without (erroneous) push_back of reference
typename boost::range_value<Container>::type p;
geometry::copy_coordinates(container.front(), p);
*(std::back_inserter(container)++) = p;
}
}
}} // namespace detail::overlay
@@ -197,146 +230,151 @@ inline void stop_gracefully(Container& container, bool& stop,
*/
template
<
typename SideStrategy,
typename GeometryOut,
typename Geometry1,
typename Geometry2,
typename TurnPoints,
typename OutputIterator
typename Turns,
typename Rings
>
inline void traverse(Geometry1 const& geometry1,
Geometry2 const& geometry2,
detail::overlay::operation_type operation,
TurnPoints& turn_points,
OutputIterator out)
Turns& turns, Rings& rings)
{
concept::check<const Geometry1>();
concept::check<const Geometry2>();
typedef typename boost::range_iterator<TurnPoints>::type turn_iterator;
typedef typename boost::range_value<TurnPoints>::type turn_type;
typedef typename boost::range_iterator<Turns>::type turn_iterator;
typedef typename boost::range_value<Turns>::type turn_type;
typedef typename boost::range_iterator
<
typename turn_type::container_type
>::type turn_operation_iterator_type;
std::size_t size_at_start = boost::size(rings);
// Iterate through all unvisited points
for (turn_iterator it = boost::begin(turn_points);
it != boost::end(turn_points);
++it)
bool fail = false;
do
{
#ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION
//std::cout << "traversal: try: " << *it;
#endif
for (turn_operation_iterator_type iit = boost::begin(it->operations);
iit != boost::end(it->operations);
++iit)
fail = false;
// Iterate through all unvisited points
for (turn_iterator it = boost::begin(turns);
! fail && it != boost::end(turns);
++it)
{
if (iit->visited.none()
&& (iit->operation == operation
|| iit->operation == detail::overlay::operation_continue)
)
// Skip the ones marked ignore (these are: "uu", so self-tangent)
if (! it->ignore)
{
set_visited_for_contine(*it, *iit);
GeometryOut current_output;
*(std::back_inserter(current_output)++) = it->point;
turn_iterator current = it;
turn_operation_iterator_type current_iit = iit;
segment_identifier current_seg_id;
#ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION
//std::cout << "traversal: start: " << *current;
#endif
detail::overlay::assign_next_ip(geometry1, geometry2,
turn_points,
current, current_output,
*iit, current_seg_id);
detail::overlay::select_next_ip(
operation,
*current,
current_seg_id,
boost::end(it->operations),
current_iit);
iit->visited.set_started();
unsigned int i = 0;
bool stop = false;
while (current_iit != iit && ! stop)
for (turn_operation_iterator_type iit = boost::begin(it->operations);
! fail && ! it->rejected && iit != boost::end(it->operations);
++iit)
{
#ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION
//std::cout << "traverse: " << *current;
#endif
if (current_iit->visited.visited())
if (iit->visited.none()
&& (iit->operation == operation
|| iit->operation == detail::overlay::operation_continue)
)
{
// It visits a visited node again, without passing the start node.
// This makes it suspicious for endless loops
// Check if it is really same node
detail::overlay::stop_gracefully(
current_output, stop, "Visit again");
}
set_visited_for_contine(*it, *iit);
typename boost::range_value<Rings>::type current_output;
*(std::back_inserter(current_output)++) = it->point;
turn_iterator current = it;
turn_operation_iterator_type current_iit = iit;
segment_identifier current_seg_id;
detail::overlay::assign_next_ip(geometry1, geometry2,
turns,
current, current_output,
*iit, current_seg_id);
if (! detail::overlay::select_next_ip(
operation,
*current,
current_seg_id,
boost::end(it->operations),
current_iit))
{
detail::overlay::backtrack(
size_at_start, fail,
rings, current_output, turns, *it,
"Dead end at start",
geometry1, geometry2);
}
else
{
iit->visited.set_started();
unsigned int i = 0;
while (current_iit != iit && ! fail)
{
if (current_iit->visited.visited())
{
// It visits a visited node again, without passing the start node.
// This makes it suspicious for endless loops
detail::overlay::backtrack(
size_at_start, fail,
rings, current_output, turns, *it,
"Visit again",
geometry1, geometry2);
}
else
{
// We assume clockwise polygons only, non self-intersecting, closed.
// However, the input might be different, and checking validity
// is up to the library user.
// We assume clockwise polygons only, non self-intersecting, closed.
// However, the input might be different, and checking validity
// is up to the library user.
// Therefore we make here some sanity checks. If the input
// violates the assumptions, the output polygon will not be correct
// but the routine will stop and output the current polygon, and
// will continue with the next one.
// Therefore we make here some sanity checks. If the input
// violates the assumptions, the output polygon will not be correct
// but the routine will stop and output the current polygon, and
// will continue with the next one.
// Below three reasons to stop.
assign_next_ip(geometry1, geometry2,
turn_points, current, current_output,
*current_iit, current_seg_id);
// Below three reasons to stop.
assign_next_ip(geometry1, geometry2,
turns, current, current_output,
*current_iit, current_seg_id);
if (! detail::overlay::select_next_ip(
operation,
*current,
current_seg_id,
iit,
current_iit))
{
// Should not occur in valid (non-self-intersecting) polygons
// Should not occur in self-intersecting polygons without spikes
// Might occur in polygons with spikes
detail::overlay::stop_gracefully(
current_output, stop, "Dead end");
}
if (! detail::overlay::select_next_ip(
operation,
*current,
current_seg_id,
iit,
current_iit))
{
// Should not occur in valid (non-self-intersecting) polygons
// Should not occur in self-intersecting polygons without spikes
// Might occur in polygons with spikes
detail::overlay::backtrack(
size_at_start, fail,
rings, current_output, turns, *it,
"Dead end",
geometry1, geometry2);
}
if (i++ > turn_points.size())
{
// Sanity check: there may be never more loops
// than intersection points.
detail::overlay::stop_gracefully(
current_output, stop, "Endless loop");
if (i++ > turns.size())
{
// Sanity check: there may be never more loops
// than intersection points.
detail::overlay::backtrack(
size_at_start, fail,
rings, current_output, turns, *it,
"Endless loop",
geometry1, geometry2);
}
}
}
if (! fail)
{
iit->visited.set_finished();
rings.push_back(current_output);
}
}
}
}
iit->visited.set_finished();
#ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION
//std::cout << "traversal: finish: " << *current;
std::cout << geometry::wkt(current_output) << std::endl;
#endif
*out++ = current_output;
}
}
}
#ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION
std::cout << "traversal: all IP's handled" << std::endl;
#endif
} while (fail);
}

2
include/boost/geometry/algorithms/union.hpp
View File

@@ -60,7 +60,7 @@ template
typename Strategy
>
struct union_inserter
: detail::overlay::overlay_and_assemble
: detail::overlay::overlay
<G1, G2, OutputIterator, GeometryOut, 1, Strategy>
{
};

4
include/boost/geometry/multi/algorithms/overlay/assemble.hpp
View File

@@ -43,7 +43,7 @@ struct add_to_containment<multi_polygon_tag, MultiPolygon>
template <typename ContainmentContainer, typename Map>
static inline void apply(ContainmentContainer& container,
ring_identifier const& id, MultiPolygon const& multi_polygon,
Map const& map)
Map const& map, bool dissolve)
{
ring_identifier copy = id;
copy.multi_index = 0;
@@ -58,7 +58,7 @@ struct add_to_containment<multi_polygon_tag, MultiPolygon>
<
polygon_tag,
typename boost::range_value<MultiPolygon>::type
>::apply(container, copy, *it, map);
>::apply(container, copy, *it, map, dissolve);
}
}
};