From c64404808cca73495659ecb486e66d94b60eb44e Mon Sep 17 00:00:00 2001 From: Menelaos Karavelas Date: Wed, 29 Apr 2015 16:20:27 +0300 Subject: [PATCH] [algorithms][envelope][spherical] move implementation for envelope(segment, mbr) to a different file and add implementation for envelope(segment, mbr) for the spherical equatorial coordinate system --- .../algorithms/detail/envelope/segment.hpp | 330 ++++++++++++++++++ 1 file changed, 330 insertions(+) create mode 100644 include/boost/geometry/algorithms/detail/envelope/segment.hpp diff --git a/include/boost/geometry/algorithms/detail/envelope/segment.hpp b/include/boost/geometry/algorithms/detail/envelope/segment.hpp new file mode 100644 index 000000000..1c546c3c5 --- /dev/null +++ b/include/boost/geometry/algorithms/detail/envelope/segment.hpp @@ -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 + +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include + + +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 + 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(); + } + + template + static inline void swap(CalculationType& lon1, + CalculationType& lat1, + CalculationType& lon2, + CalculationType& lat2) + { + std::swap(lon1, lon2); + std::swap(lat1, lat2); + } + + template + 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(); + + return (a1 < a2) + ? (a1 < pi_half && pi_half < a2) + : (a1 > pi_half && pi_half > a2); + } + + template + static inline bool crosses_antimeridian(CoordinateType const& lon1, + CoordinateType const& lon2) + { + return math::larger(math::abs(lon1 - lon2), math::pi()); + } + + template + 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 + 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 + static inline void apply(CalculationType& lon1, + CalculationType& lat1, + CalculationType& lon2, + CalculationType& lat2) + { + CalculationType const half_pi = math::half_pi(); + + 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(); + 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 + static inline void apply(Point const& p1, Point const& p2, Box& mbr) + { + typedef typename coordinate_type::type box_coordinate_type; + typedef typename fp_coordinate_type::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(lon1, lat1); + math::normalize_spheroidal_coordinates(lon2, lat2); + + compute_mbr_of_segment::apply(lon1, lat1, lon2, lat2); + + math::convert_coordinates + < + radian, typename coordinate_system::type::units + >(lon1, lat1); + + math::convert_coordinates + < + radian, typename coordinate_system::type::units + >(lon2, lat2); + + assign_values(mbr, + boost::numeric_cast(lon1), + boost::numeric_cast(lat1), + boost::numeric_cast(lon2), + boost::numeric_cast(lat2)); + } +}; + +template +struct envelope_one_segment +{ + template + static inline void apply(Point const& p1, Point const& p2, Box& mbr) + { + envelope_one_point::apply(p1, mbr); + geometry::expand(mbr, p2); + } +}; + +template <> +struct envelope_one_segment + : envelope_segment_on_spheroid +{}; + +}} // namespace detail::envelope +#endif // DOXYGEN_NO_DETAIL + + +#ifndef DOXYGEN_NO_DISPATCH +namespace dispatch +{ + + +template +struct envelope +{ + template + static inline void apply(Segment const& segment, Box& mbr) + { + typename point_type::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::type + >::apply(p[0], p[1], mbr); + } +}; + + +} // namespace dispatch +#endif + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_ENVELOPE_SEGMENT_HPP