[alsorithms] [distance] Optimization: avoid to compute vertex twice in some cases

This commit is contained in:
Vissarion Fysikopoulos 2018-04-17 12:59:52 +03:00
parent 7f5236b8fd
commit 18f1394e30
3 changed files with 74 additions and 91 deletions

View File

@ -69,6 +69,22 @@ public:
static inline bool apply(Segment const& segment,
Box const& box,
Strategy const& azimuth_strategy)
{
typedef typename point_type<Segment>::type segment_point;
segment_point vertex;
std::cout << apply(segment, box, azimuth_strategy, vertex) << std::endl;
return (apply(segment, box, azimuth_strategy, vertex) > 0);
}
// returns 0 if intersect,
// 1 if disjoint (vertex not computed),
// 2 if disjoint (vertex computed)
template <typename Segment, typename Box, typename Strategy, typename P>
static inline std::size_t apply(Segment const& segment,
Box const& box,
Strategy const& azimuth_strategy,
P& vertex)
{
assert_dimension_equal<Segment, Box>();
@ -84,7 +100,7 @@ public:
// Case 1: if box contains one of segment's endpoints then they are not disjoint
if (! disjoint_point_box(p0, box) || ! disjoint_point_box(p1, box))
{
return false;
return 0;
}
// Case 2: disjoint if bounding boxes are disjoint
@ -123,7 +139,7 @@ public:
if (disjoint_box_box(box, box_seg))
{
return true;
return 1;
}
// Case 3: test intersection by comparing angles
@ -149,7 +165,7 @@ public:
// there is an intersection
if (!(b0 && b1 && b2 && b3) && (b0 || b1 || b2 || b3))
{
return false;
return 0;
}
// Case 4: The only intersection case not covered above is when all four
@ -182,6 +198,9 @@ public:
alp1,
azimuth_strategy);
geometry::set_from_radian<0>(vertex, vertex_lon);
geometry::set_from_radian<1>(vertex, vertex_lat);
// Check if the vertex point is within the band defined by the
// minimum and maximum longitude of the box; if yes, then return
// false if the point is above the min latitude of the box; return
@ -189,11 +208,11 @@ public:
if (vertex_lon >= b_lon_min && vertex_lon <= b_lon_max
&& std::abs(vertex_lat) > std::abs(b_lat_below))
{
return false;
return 0;
}
}
return true;
return 2;
}
};

View File

@ -560,54 +560,9 @@ private:
ReturnType result;
typename other_compare<LessEqual>::type less_equal;
geometry::model::box<SegmentPoint> mbr;
//geometry::model::box<SegmentPoint> mbr;
typedef geometry::model::segment<SegmentPoint> Segment;
Segment seg(p0, p1);
geometry::envelope(seg,mbr);
typedef typename coordinate_type<SegmentPoint>::type CT;
CT lon1 = geometry::get_as_radian<0>(p0);
CT lat1 = geometry::get_as_radian<1>(p0);
CT lon2 = geometry::get_as_radian<0>(p1);
CT lat2 = geometry::get_as_radian<1>(p1);
//std::cout << "lat1=" << lat1 << " lat2=" << lat2<< std::endl;
CT vertex_lat;
CT lat_sum = lat1 + lat2;
////CT b_lat_below; //latitude of box closest to equator
if (lat_sum > CT(0))
{
vertex_lat = geometry::get_as_radian<geometry::max_corner, 1>(mbr);
//b_lat_below = b_lat_min;
} else {
vertex_lat = geometry::get_as_radian<geometry::min_corner, 1>(mbr);
//b_lat_below = b_lat_max;
}
//std::cout << "vertex lat=" << vertex_lat * geometry::math::r2d<CT>() << std::endl;
CT alp1;
geometry::strategy::azimuth::geographic<> azimuth_strategy;
azimuth_strategy.apply(lon1, lat1, lon2, lat2, alp1);
typedef typename cs_tag<Segment>::type segment_cs_type;
CT vertex_lon = geometry::formula::vertex_longitude<CT, segment_cs_type>
::apply(lon1, lat1,
lon2, lat2,
vertex_lat,
alp1,
azimuth_strategy);
//std::cout << "vertex lon=" << vertex_lon * geometry::math::r2d<CT>()
// << std::endl;
SegmentPoint p_max(vertex_lon* geometry::math::r2d<CT>(),
vertex_lat* geometry::math::r2d<CT>());
//TODO: do it for general Units
typedef geometry::model::box<BoxPoint> box;
@ -616,25 +571,63 @@ private:
geometry::get<0>(top_right),
geometry::get<1>(top_right));
//todo: do it more efficient
if (! geometry::disjoint(seg,input_box))
geometry::strategy::azimuth::geographic<> azimuth_strategy;
SegmentPoint p_max;
std::size_t disjoint_result =
geometry::detail::disjoint::
disjoint_segment_box_sphere_or_spheroid<geographic_tag>::
apply(seg,input_box,azimuth_strategy,p_max);
if (disjoint_result == 0) //intersect
{
return 0;
}
if (disjoint_result == 1) // disjoint but vertex not computed
{
typedef typename coordinate_type<SegmentPoint>::type CT;
geometry::model::box<SegmentPoint> mbr;
geometry::envelope(seg, mbr);
if (less_equal(geometry::get_as_radian<0>(bottom_left), vertex_lon))
CT lon1 = geometry::get_as_radian<0>(p0);
CT lat1 = geometry::get_as_radian<1>(p0);
CT lon2 = geometry::get_as_radian<0>(p1);
CT lat2 = geometry::get_as_radian<1>(p1);
CT vertex_lat;
CT lat_sum = lat1 + lat2;
if (lat_sum > CT(0))
{
vertex_lat = geometry::get_as_radian<geometry::max_corner, 1>(mbr);
} else {
vertex_lat = geometry::get_as_radian<geometry::min_corner, 1>(mbr);
}
CT alp1;
azimuth_strategy.apply(lon1, lat1, lon2, lat2, alp1);
typedef typename cs_tag<Segment>::type segment_cs_type;
CT vertex_lon = geometry::formula::vertex_longitude<CT, segment_cs_type>
::apply(lon1, lat1,
lon2, lat2,
vertex_lat,
alp1,
azimuth_strategy);
geometry::set_from_radian<0>(p_max, vertex_lon);
geometry::set_from_radian<1>(p_max, vertex_lat);
}
//otherwise disjoint_result == 2 i.e. disjoint and vertex computed
//inside disjoint
if (less_equal(geometry::get_as_radian<0>(bottom_left),
geometry::get_as_radian<0>(p_max)))
{
typedef cast_to_result<ReturnType> cast;
result = cast::apply(ps_strategy.apply(bottom_left, p0, p1));
}
else
{
//std::cout << "p1=" << geometry::get<0>(p1) << ","
// << geometry::get<1>(p1)
// << " p_max="
// << geometry::get<0>(p_max) << ","
// << geometry::get<1>(p_max)
// << std::endl;
result = above_of_box
<
typename other_compare<LessEqual>::type
@ -696,10 +689,6 @@ private:
{
static inline bool apply(SegmentPoint const& p0,
SegmentPoint const& p1,
BoxPoint const& bottom_left0,
BoxPoint const& top_right0,
BoxPoint const& bottom_left1,
BoxPoint const& top_right1,
BoxPoint const& corner1,
BoxPoint const& corner2,
PSStrategy const& ps_strategy,
@ -711,40 +700,19 @@ private:
>::type side;
typedef cast_to_result<ReturnType> cast;
/*
ReturnType diff0 = cast::apply(geometry::get<0>(p1))
- cast::apply(geometry::get<0>(p0));
ReturnType t_min0 = cast::apply(geometry::get<0>(bottom_left0))
- cast::apply(geometry::get<0>(p0));
ReturnType t_max0 = cast::apply(geometry::get<0>(top_right0))
- cast::apply(geometry::get<0>(p0));
*/
ReturnType diff1 = cast::apply(geometry::get<1>(p1))
- cast::apply(geometry::get<1>(p0));
// ReturnType t_min1 = cast::apply(geometry::get<1>(bottom_left1))
// - cast::apply(geometry::get<1>(p0));
// ReturnType t_max1 = cast::apply(geometry::get<1>(top_right1))
// - cast::apply(geometry::get<1>(p0));
- cast::apply(geometry::get<1>(p0));
int sign = 1;
if (diff1 < 0)
{
sign = -1;
// diff1 = -diff1;
// t_min1 = -t_min1;
// t_max1 = -t_max1;
}
// t_min0 > t_max1
//if (t_min0 * diff1 > t_max1 * diff0)
if (side::apply(p0, p1, corner1) * sign < 0)
{
result = cast::apply(ps_strategy.apply(corner1, p0, p1));
return true;
}
// t_max0 < t_min1
//if (t_max0 * diff1 < t_min1 * diff0)
if (side::apply(p0, p1, corner2) * sign > 0)
{
result = cast::apply(ps_strategy.apply(corner2, p0, p1));
@ -798,8 +766,6 @@ private:
}
if (check_generic_position::apply(p0, p1,
bottom_left, top_right,
bottom_left, top_right,
top_left, bottom_right,
ps_strategy, result))
{
@ -852,8 +818,6 @@ private:
}
if (check_generic_position::apply(p0, p1,
bottom_left, top_right,
top_right, bottom_left,
bottom_left, top_right,
ps_strategy, result))
{

View File

@ -358,6 +358,6 @@ void test_distance_segment_box(Strategy_pp const& strategy_pp,
BOOST_AUTO_TEST_CASE( test_all_point_segment )
{
test_distance_segment_box(vincenty_pp(), vincenty_ps(), vincenty_pb());
//test_distance_segment_box(thomas_pp(), thomas_ps(), thomas_pb());
test_distance_segment_box(thomas_pp(), thomas_ps(), thomas_pb());
test_distance_segment_box(andoyer_pp(), andoyer_ps(), andoyer_pb());
}