Chris@16: // Boost.Geometry (aka GGL, Generic Geometry Library)
Chris@16: 
Chris@16: // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
Chris@16: 
Chris@16: // Use, modification and distribution is subject to the Boost Software License,
Chris@16: // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
Chris@16: // http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_HUILLER_HPP
Chris@16: #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_HUILLER_HPP
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: #include <boost/geometry/strategies/spherical/distance_haversine.hpp>
Chris@16: 
Chris@16: #include <boost/geometry/core/radian_access.hpp>
Chris@16: #include <boost/geometry/util/math.hpp>
Chris@16: 
Chris@16: 
Chris@16: namespace boost { namespace geometry
Chris@16: {
Chris@16: 
Chris@16: namespace strategy { namespace area
Chris@16: {
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: /*!
Chris@16: \brief Area calculation by spherical excess / Huiller's formula
Chris@16: \ingroup strategies
Chris@16: \tparam PointOfSegment point type of segments of rings/polygons
Chris@16: \tparam CalculationType \tparam_calculation
Chris@16: \author Barend Gehrels. Adapted from:
Chris@16: - http://www.soe.ucsc.edu/~pang/160/f98/Gems/GemsIV/sph_poly.c
Chris@16: - http://williams.best.vwh.net/avform.htm
Chris@16: \note The version in Gems didn't account for polygons crossing the 180 meridian.
Chris@16: \note This version works for convex and non-convex polygons, for 180 meridian
Chris@16: crossing polygons and for polygons with holes. However, some cases (especially
Chris@16: 180 meridian cases) must still be checked.
Chris@16: \note The version which sums angles, which is often seen, doesn't handle non-convex
Chris@16: polygons correctly.
Chris@16: \note The version which sums longitudes, see
Chris@16: http://trs-new.jpl.nasa.gov/dspace/bitstream/2014/40409/1/07-03.pdf, is simple
Chris@16: and works well in most cases but not in 180 meridian crossing cases. This probably
Chris@16: could be solved.
Chris@16: 
Chris@16: \note This version is made for spherical equatorial coordinate systems
Chris@16: 
Chris@16: \qbk{
Chris@16: 
Chris@16: [heading Example]
Chris@16: [area_with_strategy]
Chris@16: [area_with_strategy_output]
Chris@16: 
Chris@16: 
Chris@16: [heading See also]
Chris@16: [link geometry.reference.algorithms.area.area_2_with_strategy area (with strategy)]
Chris@16: }
Chris@16: 
Chris@16: */
Chris@16: template
Chris@16: <
Chris@16:     typename PointOfSegment,
Chris@16:     typename CalculationType = void
Chris@16: >
Chris@16: class huiller
Chris@16: {
Chris@16: typedef typename boost::mpl::if_c
Chris@16:     <
Chris@16:         boost::is_void<CalculationType>::type::value,
Chris@16:         typename select_most_precise
Chris@16:             <
Chris@16:                 typename coordinate_type<PointOfSegment>::type,
Chris@16:                 double
Chris@16:             >::type,
Chris@16:         CalculationType
Chris@16:     >::type calculation_type;
Chris@16: 
Chris@16: protected :
Chris@16:     struct excess_sum
Chris@16:     {
Chris@16:         calculation_type sum;
Chris@16: 
Chris@16:         // Distances are calculated on unit sphere here
Chris@16:         strategy::distance::haversine<calculation_type> distance_over_unit_sphere;
Chris@16: 
Chris@16: 
Chris@16:         inline excess_sum()
Chris@16:             : sum(0)
Chris@16:             , distance_over_unit_sphere(1)
Chris@16:         {}
Chris@16:         inline calculation_type area(calculation_type radius) const
Chris@16:         {
Chris@16:             return - sum * radius * radius;
Chris@16:         }
Chris@16:     };
Chris@16: 
Chris@16: public :
Chris@16:     typedef calculation_type return_type;
Chris@16:     typedef PointOfSegment segment_point_type;
Chris@16:     typedef excess_sum state_type;
Chris@16: 
Chris@16:     inline huiller(calculation_type radius = 1.0)
Chris@16:         : m_radius(radius)
Chris@16:     {}
Chris@16: 
Chris@16:     inline void apply(PointOfSegment const& p1,
Chris@16:                 PointOfSegment const& p2,
Chris@16:                 excess_sum& state) const
Chris@16:     {
Chris@16:         if (! geometry::math::equals(get<0>(p1), get<0>(p2)))
Chris@16:         {
Chris@16:             calculation_type const half = 0.5;
Chris@16:             calculation_type const two = 2.0;
Chris@16:             calculation_type const four = 4.0;
Chris@16:             calculation_type const two_pi = two * geometry::math::pi<calculation_type>();
Chris@16:             calculation_type const half_pi = half * geometry::math::pi<calculation_type>();
Chris@16: 
Chris@16:             // Distance p1 p2
Chris@16:             calculation_type a = state.distance_over_unit_sphere.apply(p1, p2);
Chris@16: 
Chris@16:             // Sides on unit sphere to south pole
Chris@16:             calculation_type b = half_pi - geometry::get_as_radian<1>(p2);
Chris@16:             calculation_type c = half_pi - geometry::get_as_radian<1>(p1);
Chris@16: 
Chris@16:             // Semi parameter
Chris@16:             calculation_type s = half * (a + b + c);
Chris@16: 
Chris@16:             // E: spherical excess, using l'Huiller's formula
Chris@16:             // [tg(e / 4)]2   =   tg[s / 2]  tg[(s-a) / 2]  tg[(s-b) / 2]  tg[(s-c) / 2]
Chris@101:             calculation_type E = four
Chris@101:                 * atan(geometry::math::sqrt(geometry::math::abs(tan(s / two)
Chris@16:                     * tan((s - a) / two)
Chris@16:                     * tan((s - b) / two)
Chris@16:                     * tan((s - c) / two))));
Chris@16: 
Chris@16:             E = geometry::math::abs(E);
Chris@16: 
Chris@16:             // In right direction: positive, add area. In left direction: negative, subtract area.
Chris@16:             // Longitude comparisons are not so obvious. If one is negative, other is positive,
Chris@16:             // we have to take the dateline into account.
Chris@16:             // TODO: check this / enhance this, should be more robust. See also the "grow" for ll
Chris@16:             // TODO: use minmax or "smaller"/"compare" strategy for this
Chris@16:             calculation_type lon1 = geometry::get_as_radian<0>(p1) < 0
Chris@16:                 ? geometry::get_as_radian<0>(p1) + two_pi
Chris@16:                 : geometry::get_as_radian<0>(p1);
Chris@16: 
Chris@16:             calculation_type lon2 = geometry::get_as_radian<0>(p2) < 0
Chris@16:                 ? geometry::get_as_radian<0>(p2) + two_pi
Chris@16:                 : geometry::get_as_radian<0>(p2);
Chris@16: 
Chris@16:             if (lon2 < lon1)
Chris@16:             {
Chris@16:                 E = -E;
Chris@16:             }
Chris@16: 
Chris@16:             state.sum += E;
Chris@16:         }
Chris@16:     }
Chris@16: 
Chris@16:     inline return_type result(excess_sum const& state) const
Chris@16:     {
Chris@16:         return state.area(m_radius);
Chris@16:     }
Chris@16: 
Chris@16: private :
Chris@16:     /// Radius of the sphere
Chris@16:     calculation_type m_radius;
Chris@16: };
Chris@16: 
Chris@16: #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
Chris@16: 
Chris@16: namespace services
Chris@16: {
Chris@16: 
Chris@16: 
Chris@16: template <typename Point>
Chris@16: struct default_strategy<spherical_equatorial_tag, Point>
Chris@16: {
Chris@16:     typedef strategy::area::huiller<Point> type;
Chris@16: };
Chris@16: 
Chris@16: // Note: spherical polar coordinate system requires "get_as_radian_equatorial"
Chris@16: /***template <typename Point>
Chris@16: struct default_strategy<spherical_polar_tag, Point>
Chris@16: {
Chris@16:     typedef strategy::area::huiller<Point> type;
Chris@16: };***/
Chris@16: 
Chris@16: } // namespace services
Chris@16: 
Chris@16: #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
Chris@16: 
Chris@16: 
Chris@16: }} // namespace strategy::area
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: }} // namespace boost::geometry
Chris@16: 
Chris@16: #endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_HUILLER_HPP