Skip to content

Commit

Permalink
Merge branch 'bg-prepare'
Browse files Browse the repository at this point in the history
  • Loading branch information
@barendgehrels
barendgehrels committed Apr 16, 2016
2 parents f23c17b + 88d8363 commit be442ef
Show file tree
Hide file tree
Showing 11 changed files with 268 additions and 76 deletions.
1 change: 1 addition & 0 deletions doc/release_notes.qbk
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,7 @@
* [@https://svn.boost.org/trac/boost/ticket/11984 11984] union_() generates self-intersecting polygon
* [@https://svn.boost.org/trac/boost/ticket/11987 11987] rtree::remove() not compiling for geographic CS.
* [@https://svn.boost.org/trac/boost/ticket/12000 12000] Uninitialized reference in (unused) constructor of relate's mask_handler.
* [@https://svn.boost.org/trac/boost/ticket/12106 12106] Invalid assertion failure in envelope() for non-cartesian very short segments.

[*Bugfixes]

Expand Down
116 changes: 69 additions & 47 deletions include/boost/geometry/algorithms/detail/envelope/segment.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,11 @@
// Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
// Copyright (c) 2009-2015 Mateusz Loskot, London, UK.

// This file was modified by Oracle on 2015.
// Modifications copyright (c) 2015, Oracle and/or its affiliates.
// This file was modified by Oracle on 2015, 2016.
// Modifications copyright (c) 2015-2016, Oracle and/or its affiliates.

// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle

// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at
Expand Down Expand Up @@ -84,6 +85,7 @@ class compute_mbr_of_segment
private:
// computes the azimuths of the segment with endpoints (lon1, lat1)
// and (lon2, lat2)
// radians
template <typename CalculationType>
static inline void azimuths(CalculationType const& lon1,
CalculationType const& lat1,
Expand All @@ -92,7 +94,7 @@ class compute_mbr_of_segment
CalculationType& a1,
CalculationType& a2)
{
BOOST_GEOMETRY_ASSERT(math::smaller(lon1, lon2));
BOOST_GEOMETRY_ASSERT(lon1 <= lon2);

CalculationType dlon = lon2 - lon1;
CalculationType sin_dlon = sin(dlon);
Expand All @@ -108,8 +110,9 @@ class compute_mbr_of_segment
a2 = atan2(-sin_dlon * cos_lat1,
cos_lat2 * sin_lat1 - sin_lat2 * cos_lat1 * cos_dlon);
a2 += math::pi<CalculationType>();
}
}

// degrees or radians
template <typename CalculationType>
static inline void swap(CalculationType& lon1,
CalculationType& lat1,
Expand All @@ -120,6 +123,7 @@ class compute_mbr_of_segment
std::swap(lat1, lat2);
}

// radians
template <typename CalculationType>
static inline bool contains_pi_half(CalculationType const& a1,
CalculationType const& a2)
Expand All @@ -134,13 +138,20 @@ class compute_mbr_of_segment
: (a1 > pi_half && pi_half > a2);
}

template <typename CoordinateType>
// radians or degrees
template <typename Units, typename CoordinateType>
static inline bool crosses_antimeridian(CoordinateType const& lon1,
CoordinateType const& lon2)
{
return math::larger(math::abs(lon1 - lon2), math::pi<CoordinateType>());
typedef math::detail::constants_on_spheroid
<
CoordinateType, Units
> constants;

return math::abs(lon1 - lon2) > constants::half_period(); // > pi
}

// radians
template <typename CalculationType>
static inline CalculationType max_latitude(CalculationType const& azimuth,
CalculationType const& latitude)
Expand All @@ -149,74 +160,85 @@ class compute_mbr_of_segment
return acos( math::abs(cos(latitude) * sin(azimuth)) );
}

template <typename CalculationType>
// degrees or radians
template <typename Units, typename CalculationType>
static inline void compute_box_corners(CalculationType& lon1,
CalculationType& lat1,
CalculationType& lon2,
CalculationType& lat2,
CalculationType const& a1,
CalculationType const& a2)
CalculationType& lat2)
{
// coordinates are assumed to be in radians
BOOST_GEOMETRY_ASSERT(! math::larger(lon1, lon2));
BOOST_GEOMETRY_ASSERT(lon1 <= lon2);

if (math::equals(a1, a2))
CalculationType lon1_rad = math::as_radian<Units>(lon1);
CalculationType lat1_rad = math::as_radian<Units>(lat1);
CalculationType lon2_rad = math::as_radian<Units>(lon2);
CalculationType lat2_rad = math::as_radian<Units>(lat2);

CalculationType a1 = 0, a2 = 0;
azimuths(lon1_rad, lat1_rad, lon2_rad, lat2_rad, a1, a2);

if (lat1 > lat2)
{
// the segment must lie on the equator; nothing to do
BOOST_GEOMETRY_ASSERT(math::equals(lat1, CalculationType(0)));
BOOST_GEOMETRY_ASSERT(math::equals(lat2, CalculationType(0)));
return;
std::swap(lat1, lat2);
std::swap(lat1_rad, lat2_rad);
}

if (math::larger(lat1, lat2))
if (math::equals(a1, a2))
{
std::swap(lat1, lat2);
// the segment must lie on the equator or is very short
return;
}

if (contains_pi_half(a1, a2))
{
CalculationType mid_lat = lat1 + lat2;
CalculationType const mid_lat = lat1 + lat2;
if (mid_lat < 0)
{
// update using min latitude
CalculationType lat_min = -max_latitude(a1, lat1);
CalculationType const lat_min_rad = -max_latitude(a1, lat1_rad);
CalculationType const lat_min = math::from_radian<Units>(lat_min_rad);

if (math::larger(lat1, lat_min))
if (lat1 > lat_min)
{
lat1 = lat_min;
}
}
else if (mid_lat > 0)
{
// update using max latitude
CalculationType lat_max = max_latitude(a1, lat1);
CalculationType const lat_max_rad = max_latitude(a1, lat1_rad);
CalculationType const lat_max = math::from_radian<Units>(lat_max_rad);

if (math::smaller(lat2, lat_max))
if (lat2 < lat_max)
{
lat2 = lat_max;
}
}
}
}

template <typename CalculationType>
template <typename Units, typename CalculationType>
static inline void apply(CalculationType& lon1,
CalculationType& lat1,
CalculationType& lon2,
CalculationType& lat2)
{
CalculationType const half_pi = math::half_pi<CalculationType>();
typedef math::detail::constants_on_spheroid
<
CalculationType, Units
> constants;

bool is_pole1 = math::equals(math::abs(lat1), half_pi);
bool is_pole2 = math::equals(math::abs(lat2), half_pi);
bool is_pole1 = math::equals(math::abs(lat1), constants::max_latitude());
bool is_pole2 = math::equals(math::abs(lat2), constants::max_latitude());

if (is_pole1 && is_pole2)
{
// both points are poles; nothing more to do:
// longitudes are already normalized to 0
BOOST_GEOMETRY_ASSERT(lon1 == CalculationType(0)
&&
lon2 == CalculationType(0));
// but just in case
lon1 = 0;
lon2 = 0;
}
else if (is_pole1 && !is_pole2)
{
Expand All @@ -233,10 +255,10 @@ class compute_mbr_of_segment
lon2 = lon1;
}

if (math::equals(lon1, lon2))
if (lon1 == lon2)
{
// segment lies on a meridian
if (math::larger(lat1, lat2))
if (lat1 > lat2)
{
std::swap(lat1, lat2);
}
Expand All @@ -245,25 +267,22 @@ class compute_mbr_of_segment

BOOST_GEOMETRY_ASSERT(!is_pole1 && !is_pole2);

if (math::larger(lon1, lon2))
if (lon1 > lon2)
{
swap(lon1, lat1, lon2, lat2);
}

if (crosses_antimeridian(lon1, lon2))
if (crosses_antimeridian<Units>(lon1, lon2))
{
lon1 += math::two_pi<CalculationType>();
lon1 += constants::period();
swap(lon1, lat1, lon2, lat2);
}

CalculationType a1 = 0, a2 = 0;
azimuths(lon1, lat1, lon2, lat2, a1, a2);

compute_box_corners(lon1, lat1, lon2, lat2, a1, a2);
compute_box_corners<Units>(lon1, lat1, lon2, lat2);
}

public:
template <typename CalculationType, typename Box>
template <typename Units, typename CalculationType, typename Box>
static inline void apply(CalculationType lon1,
CalculationType lat1,
CalculationType lon2,
Expand All @@ -274,12 +293,12 @@ class compute_mbr_of_segment

typedef typename helper_geometry
<
Box, box_coordinate_type, radian
Box, box_coordinate_type, Units
>::type helper_box_type;

helper_box_type radian_mbr;

apply(lon1, lat1, lon2, lat2);
apply<Units>(lon1, lat1, lon2, lat2);

geometry::set
<
Expand Down Expand Up @@ -316,11 +335,14 @@ struct envelope_segment_on_sphere
Point p1_normalized = detail::return_normalized<Point>(p1);
Point p2_normalized = detail::return_normalized<Point>(p2);

compute_mbr_of_segment::apply(geometry::get_as_radian<0>(p1_normalized),
geometry::get_as_radian<1>(p1_normalized),
geometry::get_as_radian<0>(p2_normalized),
geometry::get_as_radian<1>(p2_normalized),
mbr);
typedef typename coordinate_system<Point>::type::units units_type;

compute_mbr_of_segment::template apply<units_type>(
geometry::get<0>(p1_normalized),
geometry::get<1>(p1_normalized),
geometry::get<0>(p2_normalized),
geometry::get<1>(p2_normalized),
mbr);

// now compute the envelope range for coordinates of
// dimension 2 and higher
Expand Down
17 changes: 4 additions & 13 deletions include/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -339,24 +339,15 @@ struct side_sorter
private :


std::size_t move(std::size_t index)
std::size_t move(std::size_t index) const
{
int const n = m_ranked_points.size();
int result = index + 1;
if (result >= n)
{
result = 0;
}
else if (result < 0)
{
result = n - 1;
}
return result;
std::size_t const result = index + 1;
return result >= m_ranked_points.size() ? 0 : result;
}

//! member is pointer to member (source_index or multi_index)
template <signed_size_type segment_identifier::*Member>
std::size_t move(signed_size_type member_index, std::size_t index)
std::size_t move(signed_size_type member_index, std::size_t index) const
{
std::size_t result = move(index);
while (m_ranked_points[result].seg_id.*Member != member_index)
Expand Down
61 changes: 61 additions & 0 deletions include/boost/geometry/util/math.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,8 @@
#include <boost/type_traits/is_fundamental.hpp>
#include <boost/type_traits/is_integral.hpp>

#include <boost/geometry/core/cs.hpp>

#include <boost/geometry/util/select_most_precise.hpp>

namespace boost { namespace geometry
Expand Down Expand Up @@ -601,6 +603,65 @@ inline T r2d()
}


#ifndef DOXYGEN_NO_DETAIL
namespace detail {

template <typename DegreeOrRadian>
struct as_radian
{
template <typename T>
static inline T apply(T const& value)
{
return value;
}
};

template <>
struct as_radian<degree>
{
template <typename T>
static inline T apply(T const& value)
{
return value * d2r<T>();
}
};

template <typename DegreeOrRadian>
struct from_radian
{
template <typename T>
static inline T apply(T const& value)
{
return value;
}
};

template <>
struct from_radian<degree>
{
template <typename T>
static inline T apply(T const& value)
{
return value * r2d<T>();
}
};

} // namespace detail
#endif

template <typename DegreeOrRadian, typename T>
inline T as_radian(T const& value)
{
return detail::as_radian<DegreeOrRadian>::apply(value);
}

template <typename DegreeOrRadian, typename T>
inline T from_radian(T const& value)
{
return detail::from_radian<DegreeOrRadian>::apply(value);
}


/*!
\brief Calculates the haversine of an angle
\ingroup utility
Expand Down
Loading

0 comments on commit be442ef

Please sign in to comment.