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

[algorithms][envelope][spherical] move implementation for envelope(segment, mbr)

Browse Source 
to a different file and add implementation for envelope(segment, mbr) for the
spherical equatorial coordinate system
This commit is contained in:
Menelaos Karavelas
2015-04-29 16:20:27 +03:00
parent ae08aac3a5
commit c64404808c
1 changed files with 330 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

330
include/boost/geometry/algorithms/detail/envelope/segment.hpp Normal file
View File

@@ -0,0 +1,330 @@
// Boost.Geometry (aka GGL, Generic Geometry Library)
// Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
// 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.
// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
// 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)
#ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_ENVELOPE_SEGMENT_HPP
#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_ENVELOPE_SEGMENT_HPP
#include <utility>
#include <boost/assert.hpp>
#include <boost/numeric/conversion/cast.hpp>
#include <boost/geometry/core/coordinate_system.hpp>
#include <boost/geometry/core/coordinate_type.hpp>
#include <boost/geometry/core/cs.hpp>
#include <boost/geometry/core/point_type.hpp>
#include <boost/geometry/core/radian_access.hpp>
#include <boost/geometry/core/tags.hpp>
#include <boost/geometry/util/condition.hpp>
#include <boost/geometry/util/convert_on_spheroid.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
#include <boost/geometry/algorithms/assign.hpp>
#include <boost/geometry/algorithms/expand.hpp>
#include <boost/geometry/algorithms/detail/envelope/point.hpp>
#include <boost/geometry/algorithms/dispatch/envelope.hpp>
namespace boost { namespace geometry
{
#ifndef DOXYGEN_NO_DETAIL
namespace detail { namespace envelope
{
// Computes the MBR of a segment given by (lon1, lat1) and (lon2,
// lat2), and with azimuths a1 and a2 at the two endpoints of the
// segment.
// It is assumed that the spherical coordinates of the segment are
// normalized and in radians.
// The longitudes and latitudes of the endpoints are overridden by
// those of the box.
class compute_mbr_of_segment
{
private:
// computes the azimuths of the segment with endpoints (lon1, lat1)
// and (lon2, lat2)
template <typename CalculationType>
static inline void azimuths(CalculationType const& lon1,
CalculationType const& lat1,
CalculationType const& lon2,
CalculationType const& lat2,
CalculationType& a1,
CalculationType& a2)
{
BOOST_ASSERT(math::smaller(lon1, lon2));
CalculationType dlon = lon2 - lon1;
CalculationType sin_dlon = sin(dlon);
CalculationType cos_dlon = cos(dlon);
CalculationType cos_lat1 = cos(lat1);
CalculationType cos_lat2 = cos(lat2);
CalculationType sin_lat1 = sin(lat1);
CalculationType sin_lat2 = sin(lat2);
a1 = atan2(sin_dlon * cos_lat2,
cos_lat1 * sin_lat2 - sin_lat1 * cos_lat2 * cos_dlon);
a2 = atan2(-sin_dlon * cos_lat1,
cos_lat2 * sin_lat1 - sin_lat2 * cos_lat1 * cos_dlon);
a2 += math::pi<CalculationType>();
}
template <typename CalculationType>
static inline void swap(CalculationType& lon1,
CalculationType& lat1,
CalculationType& lon2,
CalculationType& lat2)
{
std::swap(lon1, lon2);
std::swap(lat1, lat2);
}
template <typename CalculationType>
static inline bool contains_pi_half(CalculationType const& a1,
CalculationType const& a2)
{
// azimuths a1 and a2 are assumed to be in radians
BOOST_ASSERT(! math::equals(a1, a2));
static CalculationType const pi_half = math::half_pi<CalculationType>();
return (a1 < a2)
? (a1 < pi_half && pi_half < a2)
: (a1 > pi_half && pi_half > a2);
}
template <typename CoordinateType>
static inline bool crosses_antimeridian(CoordinateType const& lon1,
CoordinateType const& lon2)
{
return math::larger(math::abs(lon1 - lon2), math::pi<CoordinateType>());
}
template <typename CalculationType>
static inline CalculationType max_latitude(CalculationType const& azimuth,
CalculationType const& latitude)
{
// azimuth and latitude are assumed to be in radians
return acos( math::abs(cos(latitude) * sin(azimuth)) );
}
template <typename CalculationType>
static inline void compute_box_corners(CalculationType& lon1,
CalculationType& lat1,
CalculationType& lon2,
CalculationType& lat2,
CalculationType const& a1,
CalculationType const& a2)
{
// coordinates are assumed to be in radians
BOOST_ASSERT(! math::larger(lon1, lon2));
if (math::equals(a1, a2))
{
// the segment must lie on the equator; nothing to do
BOOST_ASSERT(math::equals(lat1, CalculationType(0)));
BOOST_ASSERT(math::equals(lat2, CalculationType(0)));
return;
}
if (math::larger(lat1, lat2))
{
std::swap(lat1, lat2);
}
if (contains_pi_half(a1, a2))
{
CalculationType mid_lat = lat1 + lat2;
if (mid_lat < 0)
{
// update using min latitude
CalculationType lat_min = -max_latitude(a1, lat1);
if (math::larger(lat1, lat_min))
{
lat1 = lat_min;
}
}
else if (mid_lat > 0)
{
// update using max latitude
CalculationType lat_max = max_latitude(a1, lat1);
if (math::smaller(lat2, lat_max))
{
lat2 = lat_max;
}
}
}
}
public:
template <typename CalculationType>
static inline void apply(CalculationType& lon1,
CalculationType& lat1,
CalculationType& lon2,
CalculationType& lat2)
{
CalculationType const half_pi = math::half_pi<CalculationType>();
bool is_pole1 = math::equals(math::abs(lat1), half_pi);
bool is_pole2 = math::equals(math::abs(lat2), half_pi);
if (is_pole1 && is_pole2)
{
// both points are poles; nothing more to do:
// longitudes are already normalized to 0
BOOST_ASSERT(lon1 == CalculationType(0)
&&
lon2 == CalculationType(0));
}
else if (is_pole1 && !is_pole2)
{
// first point is a pole, second point is not:
// make the longitude of the first point the same as that
// of the second point
lon1 = lon2;
}
else if (!is_pole1 && is_pole2)
{
// second point is a pole, first point is not:
// make the longitude of the second point the same as that
// of the first point
lon2 = lon1;
}
if (math::equals(lon1, lon2))
{
// segment lies on a meridian
if (math::larger(lat1, lat2))
{
std::swap(lat1, lat2);
}
return;
}
BOOST_ASSERT(!is_pole1 && !is_pole2);
if (math::larger(lon1, lon2))
{
swap(lon1, lat1, lon2, lat2);
}
if (crosses_antimeridian(lon1, lon2))
{
lon1 += math::two_pi<CalculationType>();
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);
}
};
struct envelope_segment_on_spheroid
{
template <typename Point, typename Box>
static inline void apply(Point const& p1, Point const& p2, Box& mbr)
{
typedef typename coordinate_type<Box>::type box_coordinate_type;
typedef typename fp_coordinate_type<Point>::type calculation_type;
calculation_type lon1 = geometry::get_as_radian<0>(p1);
calculation_type lat1 = geometry::get_as_radian<1>(p1);
calculation_type lon2 = geometry::get_as_radian<0>(p2);
calculation_type lat2 = geometry::get_as_radian<1>(p2);
math::normalize_spheroidal_coordinates<radian>(lon1, lat1);
math::normalize_spheroidal_coordinates<radian>(lon2, lat2);
compute_mbr_of_segment::apply(lon1, lat1, lon2, lat2);
math::convert_coordinates
<
radian, typename coordinate_system<Box>::type::units
>(lon1, lat1);
math::convert_coordinates
<
radian, typename coordinate_system<Box>::type::units
>(lon2, lat2);
assign_values(mbr,
boost::numeric_cast<box_coordinate_type>(lon1),
boost::numeric_cast<box_coordinate_type>(lat1),
boost::numeric_cast<box_coordinate_type>(lon2),
boost::numeric_cast<box_coordinate_type>(lat2));
}
};
template <typename CSTag>
struct envelope_one_segment
{
template<typename Point, typename Box>
static inline void apply(Point const& p1, Point const& p2, Box& mbr)
{
envelope_one_point<CSTag>::apply(p1, mbr);
geometry::expand(mbr, p2);
}
};
template <>
struct envelope_one_segment<spherical_equatorial_tag>
: envelope_segment_on_spheroid
{};
}} // namespace detail::envelope
#endif // DOXYGEN_NO_DETAIL
#ifndef DOXYGEN_NO_DISPATCH
namespace dispatch
{
template <typename Segment>
struct envelope<Segment, segment_tag>
{
template <typename Box>
static inline void apply(Segment const& segment, Box& mbr)
{
typename point_type<Segment>::type p[2];
detail::assign_point_from_index<0>(segment, p[0]);
detail::assign_point_from_index<1>(segment, p[1]);
detail::envelope::envelope_one_segment
<
typename cs_tag<Segment>::type
>::apply(p[0], p[1], mbr);
}
};
} // namespace dispatch
#endif
}} // namespace boost::geometry
#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_ENVELOPE_SEGMENT_HPP