boost_mirror/geometry
1
0
Fork 0
mirror of https://github.com/boostorg/geometry.git synced 2025-05-09 15:14:02 +00:00
Code Issues Packages Projects Releases Wiki Activity

[fix] for integer coordinates, segment intersection now rounds to nearest point

Browse Source
This commit is contained in:
Barend Gehrels 2022-01-26 13:48:23 +01:00
parent 6099339544
commit 9c18680f31
12 changed files with 169 additions and 75 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
include/boost/geometry/algorithms/detail/overlay/overlay.hpp
View File

@ -369,12 +369,9 @@ std::cout << "traverse" << std::endl;
// Add rings created during traversal
{
ring_identifier id(2, 0, -1);
for (typename boost::range_iterator<ring_container_type>::type
it = boost::begin(rings);
it != boost::end(rings);
++it)
for (auto const& ring : rings)
{
selected_ring_properties[id] = properties(*it, strategy);
selected_ring_properties[id] = properties(ring, strategy);
selected_ring_properties[id].reversed = ReverseOut;
id.multi_index++;
}

4
include/boost/geometry/core/config.hpp
View File

@ -35,8 +35,8 @@
// Rescaling is turned on, unless NO_ROBUSTNESS is defined
// In future versions of Boost.Geometry, it will be turned off by default
#if ! defined(BOOST_GEOMETRY_NO_ROBUSTNESS)
#define BOOST_GEOMETRY_USE_RESCALING
#if !defined(BOOST_GEOMETRY_NO_ROBUSTNESS) && !defined(BOOST_GEOMETRY_USE_RESCALING)
#define BOOST_GEOMETRY_USE_RESCALING
#endif
#endif // BOOST_GEOMETRY_CORE_CONFIG_HPP

13
include/boost/geometry/strategies/cartesian/intersection.hpp
View File

@ -170,9 +170,8 @@ struct cartesian_segments
SegmentRatio const& ratio) const
{
// Calculate the intersection point based on segment_ratio
// Up to now, division was postponed. Here we divide using numerator/
// denominator. In case of integer this results in an integer
// division.
// The division, postponed until here, is done now. In case of integer this
// results in an integer which rounds to the nearest integer.
BOOST_GEOMETRY_ASSERT(ratio.denominator() != typename SegmentRatio::int_type(0));
typedef typename promote_integral<CoordinateType>::type calc_type;
@ -185,11 +184,11 @@ struct cartesian_segments
calc_type const dy_calc = boost::numeric_cast<calc_type>(dy);
set<0>(point, get<0, 0>(segment)
+ boost::numeric_cast<CoordinateType>(numerator * dx_calc
/ denominator));
+ boost::numeric_cast<CoordinateType>(
math::divide<calc_type>(numerator * dx_calc, denominator)));
set<1>(point, get<0, 1>(segment)
+ boost::numeric_cast<CoordinateType>(numerator * dy_calc
/ denominator));
+ boost::numeric_cast<CoordinateType>(
math::divide<calc_type>(numerator * dy_calc, denominator)));
}
template <int Index, int Dim, typename Point, typename Segment>

35
include/boost/geometry/util/math.hpp
View File

@ -525,6 +525,30 @@ struct rounding_cast<Result, Source, true, false>
}
};
template <typename T, bool IsIntegral = std::is_integral<T>::value>
struct divide
{
static inline T apply(T const& n, T const& d)
{
return n / d;
}
};
template <typename T>
struct divide<T, true>
{
static inline T apply(T const& n, T const& d)
{
return n == 0 ? 0
: n < 0
? (d < 0 ? (n + (-d + 1) / 2) / d + 1
: (n + ( d + 1) / 2) / d - 1 )
: (d < 0 ? (n - (-d + 1) / 2) / d - 1
: (n - ( d + 1) / 2) / d + 1 )
;
}
};
} // namespace detail
#endif
@ -796,6 +820,17 @@ inline Result rounding_cast(T const& v)
return detail::rounding_cast<Result, T>::apply(v);
}
/*
\brief Short utility to divide. If the division is integer, it rounds the division
to the nearest value, without using floating point calculations
\ingroup utility
*/
template <typename T>
inline T divide(T const& n, T const& d)
{
return detail::divide<T>::apply(n, d);
}
/*!
\brief Evaluate the sine and cosine function with the argument in degrees
\note The results obey exactly the elementary properties of the trigonometric

7
test/algorithms/overlay/multi_overlay_cases.hpp
View File

File diff suppressed because one or more lines are too long

55
test/algorithms/overlay/overlay.cpp
View File

@ -18,10 +18,6 @@
#include <sstream>
#include <string>
#if defined(TEST_WITH_SVG)
# include <boost/geometry/io/svg/svg_mapper.hpp>
#endif
#include <geometry_test_common.hpp>
#include <algorithms/check_validity.hpp>
@ -29,15 +25,12 @@
#include <boost/geometry/algorithms/detail/overlay/debug_turn_info.hpp>
#include <boost/geometry/algorithms/detail/overlay/overlay.hpp>
#include <boost/geometry/geometries/geometries.hpp>
//#include <boost/geometry/extensions/algorithms/inverse.hpp>
#include <boost/geometry/io/wkt/wkt.hpp>
#if defined(TEST_WITH_SVG)
# include <boost/geometry/io/svg/svg_mapper.hpp>
#endif
#include <boost/geometry/io/wkt/read.hpp>
#include "multi_overlay_cases.hpp"
@ -71,9 +64,8 @@ struct map_visitor
template <typename Turns>
void visit_turns(int phase, Turns const& turns)
{
typedef typename boost::range_value<Turns>::type turn_type;
int index = 0;
for (turn_type const& turn : turns)
for (auto const& turn : turns)
{
switch (phase)
{
@ -128,9 +120,8 @@ struct map_visitor
template <typename Clusters, typename Turns>
void visit_clusters(Clusters const& clusters, Turns const& turns)
{
typedef typename boost::range_value<Turns>::type turn_type;
int index = 0;
for (turn_type const& turn : turns)
for (auto const& turn : turns)
{
if (turn.cluster_id >= 0)
{
@ -158,8 +149,6 @@ struct map_visitor
template <typename Turns, typename Turn, typename Operation>
void visit_traverse(Turns const& turns, Turn const& turn, Operation const& op, const std::string& header)
{
typedef typename boost::range_value<Turns>::type turn_type;
if (! m_do_output)
{
return;
@ -381,7 +370,7 @@ void test_overlay(std::string const& caseid,
std::ofstream svg(filename.str().c_str());
typedef bg::svg_mapper<typename bg::point_type<Geometry>::type> svg_mapper;
using svg_mapper = bg::svg_mapper<typename bg::point_type<Geometry>::type>;
svg_mapper mapper(svg, 500, 500);
mapper.add(g1);
@ -395,8 +384,8 @@ void test_overlay(std::string const& caseid,
#endif
typedef typename boost::range_value<Geometry>::type geometry_out;
typedef bg::detail::overlay::overlay
using geometry_out = typename boost::range_value<Geometry>::type ;
using overlay = bg::detail::overlay::overlay
<
Geometry, Geometry,
bg::detail::overlay::do_reverse<bg::point_order<Geometry>::value>::value,
@ -406,20 +395,20 @@ void test_overlay(std::string const& caseid,
bg::detail::overlay::do_reverse<bg::point_order<Geometry>::value>::value,
geometry_out,
OverlayType
> overlay;
>;
typedef typename bg::strategies::relate::services::default_strategy
using strategy_type = typename bg::strategies::relate::services::default_strategy
<
Geometry, Geometry
>::type strategy_type;
>::type;
strategy_type strategy;
typedef typename bg::rescale_overlay_policy_type
using rescale_policy_type = typename bg::rescale_overlay_policy_type
<
Geometry,
Geometry
>::type rescale_policy_type;
>::type;
rescale_policy_type robust_policy
= bg::get_rescale_policy<rescale_policy_type>(g1, g2);
@ -474,9 +463,9 @@ void test_overlay(std::string const& caseid,
template <typename T, bool Clockwise>
void test_all()
{
typedef bg::model::point<T, 2, bg::cs::cartesian> point_type;
typedef bg::model::polygon<point_type, Clockwise> polygon;
typedef bg::model::multi_polygon<polygon> multi_polygon;
using point_type = bg::model::point<T, 2, bg::cs::cartesian>;
using polygon = bg::model::polygon<point_type, Clockwise>;
using multi_polygon = bg::model::multi_polygon<polygon>;
TEST_UNION(case_multi_simplex, 14.58, 1, 0);
TEST_INTERSECTION(case_multi_simplex, 6.42, 2, 0);
@ -510,16 +499,22 @@ void test_all()
TEST_UNION(case_recursive_boxes_12, 6.0, 6, 0);
TEST_UNION(case_recursive_boxes_13, 10.25, 3, 0);
TEST_INTERSECTION(issue_930, 8.3333333, 1, 0);
}
// std::cout
// << " \""
// << bg::inverse<multi_polygon>(case_65_multi[0], 1.0)
// << "\"" << std::endl;
template <typename T, bool Clockwise>
void test_integer()
{
using point_type = bg::model::point<T, 2, bg::cs::cartesian>;
using polygon = bg::model::polygon<point_type, Clockwise>;
using multi_polygon = bg::model::multi_polygon<polygon>;
TEST_INTERSECTION(issue_930, 10, 1, 0);
}
int test_main(int, char* [])
{
test_integer<int, true>();
test_all<double, true>();
// test_all<double, false>();
return 0;
}

13
test/algorithms/set_operations/difference/difference.cpp
View File

@ -612,9 +612,9 @@ void test_specific()
test_one<polygon, polygon, polygon>("ggl_list_20120717_volker",
ggl_list_20120717_volker[0], ggl_list_20120717_volker[1],
1, 11, 3371540,
1, 4, 385,
1, 16, 3371540 + 385);
1, 11, 3370444,
1, 4, 383,
2, 16, 3370444 + 383);
test_one<polygon, polygon, polygon>("ticket_10658",
ticket_10658[0], ticket_10658[1],
@ -623,11 +623,12 @@ void test_specific()
test_one<polygon, polygon, polygon>("ticket_11121",
ticket_11121[0], ticket_11121[1],
2, 8, 489763.5,
1, 4, 6731652.0);
2, 8, 489904.5,
1, 4, 6755355,
1);
// Generates spikes, both a-b and b-a
TEST_DIFFERENCE(ticket_11676, 2, 2537992.5, 2, 294963.5, 3);
TEST_DIFFERENCE(ticket_11676, 2, 2537404.5, 2, 295353, 3);
}
int test_main(int, char* [])

4
test/algorithms/set_operations/difference/difference_areal_linear.cpp
View File

@ -240,8 +240,8 @@ int test_main(int, char* [])
test_all<bg::model::d2::point_xy<double> >();
test_ticket_10835<int>
("MULTILINESTRING((5239 2113,5232 2115),(4794 2205,1020 2986))",
"MULTILINESTRING((5239 2113,5232 2115),(4794 2205,1460 2895))");
("MULTILINESTRING((5239 2113,5233 2114),(4794 2205,1020 2986))",
"MULTILINESTRING((5239 2113,5233 2114),(4794 2205,1460 2895))");
test_ticket_10835<double>
("MULTILINESTRING((5239 2113,5232.52 2114.34),(4794.39 2205,1020 2986))",

27
test/algorithms/set_operations/difference/difference_multi.cpp
View File

@ -440,10 +440,10 @@ void test_specific_areal()
settings.set_test_validity(false);
TEST_DIFFERENCE_WITH(0, 1, ticket_11674,
count_set(3, 4),
expectation_limits(9105473, 9105782),
5,
expectation_limits(119059, 119424), ignore_count());
count_set(2, 3, 4),
expectation_limits(9105196, 9105705),
6,
expectation_limits(119059, 119757), ignore_count());
}
{
@ -458,8 +458,8 @@ void test_specific_areal()
expectation_limits(597.0, 598.0), 2);
TEST_DIFFERENCE_WITH(2, 3, ticket_12751,
2, expectation_limits(2537992, 2538306),
2, expectation_limits(294736, 294964),
2, expectation_limits(2537404, 2538306),
2, expectation_limits(294736, 295353),
3);
}
@ -469,9 +469,10 @@ void test_specific_areal()
ut_settings settings;
settings.remove_spikes = true;
settings.sym_difference = false;
settings.set_test_validity(false);
TEST_DIFFERENCE_WITH(0, 1, ticket_12752,
count_set(2, 3), expectation_limits(2776656, 2776693),
3, expectation_limits(7710, 7894),
count_set(2, 3), expectation_limits(2775740, 2776693),
3, expectation_limits(7710, 7903),
2);
}
@ -479,13 +480,13 @@ void test_specific_areal()
const std::string a_min_b =
test_one<Polygon, MultiPolygon, MultiPolygon>("ticket_10661",
ticket_10661[0], ticket_10661[1],
2, -1, expectation_limits(1441632, 1441855),
2, -1, expectation_limits(1441855, 1442081),
2, -1, expectation_limits(13167454, 13182415),
count_set(3, 4));
test_one<Polygon, MultiPolygon, MultiPolygon>("ticket_10661_2",
a_min_b, ticket_10661[2],
1, 8, 825192.0,
1, 8, 825640.5,
1, 10, expectation_limits(27226148, 27842812),
count_set(1, 2));
}
@ -495,9 +496,9 @@ void test_specific_areal()
settings.sym_difference = false;
TEST_DIFFERENCE_WITH(0, 1, ticket_9942, 4,
expectation_limits(7427727, 7428108), 4,
expectation_limits(130083, 131823), 4);
TEST_DIFFERENCE_WITH(0, 1, ticket_9942a, 2,
expectation_limits(7427727, 7428549), 5,
expectation_limits(184653, 184687), 5);
TEST_DIFFERENCE_WITH(0, 1, ticket_9942a, 1,
expectation_limits(412676, 413469), 2,
expectation_limits(76779, 77038), 4);
}

16
test/algorithms/set_operations/difference/difference_multi_spike.cpp
View File

@ -41,17 +41,17 @@ void test_spikes_in_ticket_8364()
test_one<polygon, multi_polygon, multi_polygon>("ticket_8364_step3",
"MULTIPOLYGON(((3232 2532,2136 2790,1032 1764,1032 1458,1032 1212,2136 2328,3232 2220,3232 1056,1031 1056,1031 2856,3232 2856,3232 2532)))",
"MULTIPOLYGON(((1032 2130,2052 2712,1032 1764,1032 2130)),((3234 2580,3234 2532,2558 2690,3234 2580)),((2558 2690,2136 2760,2052 2712,2136 2790,2558 2690)))",
count_set(2, 3), -1, expectation_limits(2775256, 2775610), // SQL Server: 2775256.47588724
3, -1, expectation_limits(7810, 7893), // SQL Server: 7810.48711165739
count_set(2, 6), ignore_validity);
count_set(1, 2, 3), -1, expectation_limits(2774644, 2775610), // SQL Server: 2775256.47588724
3, -1, expectation_limits(7810, 7903), // SQL Server: 7810.48711165739
count_set(2, 3, 6), ignore_validity);
// TODO: behaviour is not good yet. It is changed at introduction of self-turns.
test_one<polygon, multi_polygon, multi_polygon>("ticket_8364_step4",
"MULTIPOLYGON(((2567 2688,2136 2790,2052 2712,1032 2130,1032 1764,1032 1458,1032 1212,2136 2328,3232 2220,3232 1056,1031 1056,1031 2856,3232 2856,3232 2580,2567 2688)))",
"MULTIPOLYGON(((1032 2556,1778 2556,1032 2130,1032 2556)),((3234 2580,3234 2556,1778 2556,2136 2760,3234 2580)))",
count_set(1, 2), -1, expectation_limits(2615783, 2616030), // SQL Server: 2616029.55616044
1, -1, expectation_limits(161054, 161134), // SQL Server: 161054.560110092
count_set(1, 3), ignore_validity);
count_set(1, 2), -1, expectation_limits(2615783, 2616400), // SQL Server: 2616029.55616044
1, -1, expectation_limits(160954, 161134), // SQL Server: 161054.560110092
count_set(1, 2, 3), ignore_validity);
}
template <typename P, bool ClockWise, bool Closed>
@ -75,10 +75,6 @@ void test_spikes_in_ticket_8365()
count_set(2, 4));
}
int test_main(int, char* [])
{
test_spikes_in_ticket_8364<bg::model::d2::point_xy<double>, true, true>();

1
test/util/Jamfile
View File

@ -20,6 +20,7 @@ test-suite boost-geometry-util
[ run calculation_type.cpp : : : : util_calculation_type ]
[ run for_each_coordinate.cpp : : : : util_for_each_coordinate ]
[ run math_abs.cpp : : : : util_math_abs ]
[ run math_divide.cpp : : : : util_math_divide ]
[ run math_equals.cpp : : : : util_math_equals ]
[ run math_sqrt.cpp : : : : util_math_sqrt ]
[ run promote_integral.cpp : : : : util_promote_integral ]

62
test/util/math_divide.cpp Normal file
View File

@ -0,0 +1,62 @@
// Boost.Geometry
// Unit Test
// Copyright (c) 2022 Barend Gehrels, 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)
#include <geometry_test_common.hpp>
#include <boost/geometry/util/math.hpp>
namespace bgm = bg::math;
template <typename T>
void test_integer()
{
BOOST_CHECK_EQUAL(0, bgm::divide<T>(0, 2));
BOOST_CHECK_EQUAL(1, bgm::divide<T>(1, 2));
BOOST_CHECK_EQUAL(-1, bgm::divide<T>(-1, 2));
BOOST_CHECK_EQUAL(-1, bgm::divide<T>(1, -2));
BOOST_CHECK_EQUAL(1, bgm::divide<T>(-1, -2));
BOOST_CHECK_EQUAL(2, bgm::divide<T>(4, 2));
BOOST_CHECK_EQUAL(3, bgm::divide<T>(5, 2));
BOOST_CHECK_EQUAL(9, bgm::divide<T>(94, 10));
BOOST_CHECK_EQUAL(10, bgm::divide<T>(95, 10));
BOOST_CHECK_EQUAL(10, bgm::divide<T>(99, 10));
BOOST_CHECK_EQUAL(-9, bgm::divide<T>(94, -10));
BOOST_CHECK_EQUAL(-10, bgm::divide<T>(95, -10));
BOOST_CHECK_EQUAL(-10, bgm::divide<T>(99, -10));
BOOST_CHECK_EQUAL(-9, bgm::divide<T>(-94, 10));
BOOST_CHECK_EQUAL(-10, bgm::divide<T>(-95, 10));
BOOST_CHECK_EQUAL(-10, bgm::divide<T>(-99, 10));
BOOST_CHECK_EQUAL(9, bgm::divide<T>(-94, -10));
BOOST_CHECK_EQUAL(10, bgm::divide<T>(-95, -10));
BOOST_CHECK_EQUAL(10, bgm::divide<T>(-99, -10));
BOOST_CHECK_EQUAL(5, bgm::divide<T>(4567, 1000));
}
template <typename T>
void test_floating_point()
{
constexpr double eps = 1.0e-6;
BOOST_CHECK_CLOSE(4.567, bgm::divide<T>(4567, 1000), eps);
}
int test_main(int, char* [])
{
test_integer<int>();
test_floating_point<double>();
return 0;
}