[Boost-commit] svn:boost r63839 - in sandbox/geometry/boost/geometry: extensions/contrib extensions/gis/geographic/strategies util

Subject: [Boost-commit] svn:boost r63839 - in sandbox/geometry/boost/geometry: extensions/contrib extensions/gis/geographic/strategies util
From: barend.gehrels_at_[hidden]
Date: 2010-07-11 05:40:33

Author: barendgehrels
Date: 2010-07-11 05:40:32 EDT (Sun, 11 Jul 2010)
New Revision: 63839
URL: http://svn.boost.org/trac/boost/changeset/63839

Log:
Updated definition of PI to support templated UDT 
Updated Andoyer for high precision
Text files modified: 
   sandbox/geometry/boost/geometry/extensions/contrib/ttmath_stub.hpp                |    45 ++++++------                            
   sandbox/geometry/boost/geometry/extensions/gis/geographic/strategies/andoyer.hpp  |    63 ++++++++++++------                      
   sandbox/geometry/boost/geometry/extensions/gis/geographic/strategies/vincenty.hpp |   133 +++++++++++++++++++-------------------- 
   sandbox/geometry/boost/geometry/util/math.hpp                                     |    30 +++++++-                                
   4 files changed, 151 insertions(+), 120 deletions(-)

Modified: sandbox/geometry/boost/geometry/extensions/contrib/ttmath_stub.hpp
==============================================================================
--- sandbox/geometry/boost/geometry/extensions/contrib/ttmath_stub.hpp	(original)
+++ sandbox/geometry/boost/geometry/extensions/contrib/ttmath_stub.hpp	2010-07-11 05:40:32 EDT (Sun, 11 Jul 2010)
@@ -2,6 +2,7 @@
 #define TTMATH_STUB
 
 #include <boost/math/constants/constants.hpp>
+#include <boost/geometry/util/math.hpp>
 
 
 #include <ttmath/ttmath.h>
@@ -107,39 +108,39 @@
     */
 };
 
-namespace boost{ namespace math { namespace constants
+namespace boost{ namespace geometry { namespace math
+{
+
+namespace detail
 {
     // Workaround for boost::math::constants::pi:
     // 1) lexical cast -> stack overflow and
     // 2) because it is implemented as a function, generic implementation not possible
 
+    // Partial specialization for ttmath
     template <ttmath::uint Exponent, ttmath::uint Mantissa>
-    inline ttmath::Big<Exponent, Mantissa> ttmath_pi()
+    struct define_pi<ttmath::Big<Exponent, Mantissa> >
     {
-        static ttmath::Big<Exponent, Mantissa> const the_pi("3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196");
-        return the_pi;
-    }
-
+        static inline ttmath::Big<Exponent, Mantissa> apply()
+        {
+            static ttmath::Big<Exponent, Mantissa> const the_pi(
+                "3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196");
+            return the_pi;
+        }
+    };
 
-    template <> inline ttmath::Big<1,4> pi<ttmath::Big<1,4> >(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC( (ttmath::Big<1,4>) ))
-    {
-        return ttmath_pi<1,4>();
-    }
-    template <> inline ttmath::Big<1,8> pi<ttmath::Big<1,8> >(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC( (ttmath::Big<1,8>) ))
+    template <>
+    struct define_pi<ttmath_big>
     {
-        return ttmath_pi<1,8>();
-    }
-
-    // and so on...
-    // or fix lexical_cast but probably has the same problem
-
+        static inline ttmath_big apply()
+        {
+            return define_pi<ttmath::Big<1,4> >::apply();
+        }
+    };
 
-    template <> inline ttmath_big pi<ttmath_big >(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(ttmath_big))
-    {
-        return ttmath_pi<1,4>();
-    }
+}
+}}} // boost::geometry::math
 
-}}}
 
 
 #endif

Modified: sandbox/geometry/boost/geometry/extensions/gis/geographic/strategies/andoyer.hpp
==============================================================================
--- sandbox/geometry/boost/geometry/extensions/gis/geographic/strategies/andoyer.hpp	(original)
+++ sandbox/geometry/boost/geometry/extensions/gis/geographic/strategies/andoyer.hpp	2010-07-11 05:40:32 EDT (Sun, 11 Jul 2010)
@@ -13,6 +13,9 @@
 #include <boost/geometry/strategies/distance.hpp>
 #include <boost/geometry/core/radian_access.hpp>
 #include <boost/geometry/core/coordinate_type.hpp>
+#include <boost/geometry/util/select_calculation_type.hpp>
+#include <boost/geometry/util/promote_floating_point.hpp>
+#include <boost/geometry/util/math.hpp>
 
 #include <boost/geometry/extensions/gis/geographic/detail/ellipsoid.hpp>
 
@@ -27,8 +30,8 @@
 /*!
     \brief Point-point distance approximation taking flattening into account
     \ingroup distance
-    \tparam P1 first point type
-    \tparam P2 optional second point type
+    \tparam Point1 first point type
+    \tparam Point2 optional second point type
     \author After Andoyer, 19xx, republished 1950, republished by Meeus, 1999
     \note Although not so well-known, the approximation is very good: in all cases the results
     are about the same as Vincenty. In my (Barend's) testcases the results didn't differ more than 6 m
@@ -41,20 +44,28 @@
 */
 template
 <
-    typename P1,
-    typename P2 = P1
-    // calculation_type
+    typename Point1,
+    typename Point2 = Point1,
+    typename CalculationType = void
 >
 class andoyer
 {
     public :
-        typedef double return_type;
+    typedef typename promote_floating_point
+        <
+            typename select_calculation_type
+                <
+                    Point1,
+                    Point2,
+                    CalculationType
+                >::type
+        >::type calculation_type;
 
-        andoyer()
+        inline andoyer()
             : m_ellipsoid()
         {}
 
-        explicit inline andoyer(double f)
+        explicit inline andoyer(calculation_type f)
             : m_ellipsoid(f)
         {}
 
@@ -63,10 +74,10 @@
         {}
 
 
-        inline return_type apply(P1 const& p1, P2 const& p2) const
+        inline calculation_type apply(Point1 const& point1, Point2 const& point2) const
         {
-            return calc(get_as_radian<0>(p1), get_as_radian<1>(p1),
-                            get_as_radian<0>(p2), get_as_radian<1>(p2));
+            return calc(get_as_radian<0>(point1), get_as_radian<1>(point1),
+                            get_as_radian<0>(point2), get_as_radian<1>(point2));
         }
 
         inline geometry::detail::ellipsoid ellipsoid() const
@@ -76,13 +87,15 @@
 
 
     private :
-        typedef typename coordinate_type<P1>::type T1;
-        typedef typename coordinate_type<P2>::type T2;
+        typedef typename coordinate_type<Point1>::type T1;
+        typedef typename coordinate_type<Point2>::type T2;
         geometry::detail::ellipsoid m_ellipsoid;
 
-        inline return_type calc(T1 const& lon1, T1 const& lat1, T2 const& lon2, T2 const& lat2) const
+        inline calculation_type calc(calculation_type const& lon1, 
+                    calculation_type const& lat1, 
+                    calculation_type const& lon2, 
+                    calculation_type const& lat2) const
         {
-            typedef double calculation_type;
             calculation_type G = (lat1 - lat2) / 2.0;
             calculation_type lambda = (lon1 - lon2) / 2.0;
 
@@ -109,16 +122,20 @@
                 return 0.0;
             }
 
+            calculation_type const c1 = 1;
+            calculation_type const c2 = 2;
+            calculation_type const c3 = 3;
+
             calculation_type omega = atan(sqrt(S / C));
             calculation_type r = sqrt(S * C) / omega; // not sure if this is r or greek nu
 
-            calculation_type D = 2.0 * omega * m_ellipsoid.a();
-            calculation_type H1 = (3 * r - 1.0) / (2.0 * C);
-            calculation_type H2 = (3 * r + 1.0) / (2.0 * S);
-
-            return return_type(D
-                * (1.0 + m_ellipsoid.f() * H1 * sinF2 * cosG2
-                            - m_ellipsoid.f() * H2 * cosF2 * sinG2));
+            calculation_type D = c2 * omega * m_ellipsoid.a();
+            calculation_type H1 = (c3 * r - c1) / (c2 * C);
+            calculation_type H2 = (c3 * r + c1) / (c2 * S);
+
+            calculation_type f = m_ellipsoid.f();
+
+            return D * (c1 + f * H1 * sinF2 * cosG2 - f * H2 * cosF2 * sinG2);
         }
 };
 
@@ -137,7 +154,7 @@
 template <typename Point1, typename Point2>
 struct return_type<strategy::distance::andoyer<Point1, Point2> >
 {
-    typedef typename strategy::distance::andoyer<Point1, Point2>::return_type type;
+    typedef typename strategy::distance::andoyer<Point1, Point2>::calculation_type type;
 };
 
 

Modified: sandbox/geometry/boost/geometry/extensions/gis/geographic/strategies/vincenty.hpp
==============================================================================
--- sandbox/geometry/boost/geometry/extensions/gis/geographic/strategies/vincenty.hpp	(original)
+++ sandbox/geometry/boost/geometry/extensions/gis/geographic/strategies/vincenty.hpp	2010-07-11 05:40:32 EDT (Sun, 11 Jul 2010)
@@ -61,9 +61,9 @@
                             Point2,
                             CalculationType
                         >::type,
-                    double // it should be at least double otherwise Vincenty does not run
+                    double // to avoid bad results in float
                 >::type
-        >::type return_type;
+        >::type calculation_type;
 
     inline vincenty()
     {}
@@ -72,7 +72,7 @@
         : m_ellipsoid(e)
     {}
 
-    inline return_type apply(Point1 const& p1, Point2 const& p2) const
+    inline calculation_type apply(Point1 const& p1, Point2 const& p2) const
     {
         return calculate(get_as_radian<0>(p1), get_as_radian<1>(p1),
                         get_as_radian<0>(p2), get_as_radian<1>(p2));
@@ -85,111 +85,106 @@
 
 
 private :
-    typedef return_type promoted_type;
     geometry::detail::ellipsoid m_ellipsoid;
 
-    inline return_type calculate(promoted_type const& lon1,
-                promoted_type const& lat1,
-                promoted_type const& lon2,
-                promoted_type const& lat2) const
-    {
-        namespace mc = boost::math::constants;
-
-        promoted_type const c2 = 2;
-        promoted_type const two_pi = c2 * mc::pi<promoted_type>();
+    inline calculation_type calculate(calculation_type const& lon1,
+                calculation_type const& lat1,
+                calculation_type const& lon2,
+                calculation_type const& lat2) const
+    {
+        calculation_type const c2 = 2;
+        calculation_type const pi = geometry::math::pi<calculation_type>();
+        calculation_type const two_pi = c2 * pi;
 
         // lambda: difference in longitude on an auxiliary sphere
-        promoted_type L = lon2 - lon1;
-        promoted_type lambda = L;
+        calculation_type L = lon2 - lon1;
+        calculation_type lambda = L;
 
-        if (L < -mc::pi<promoted_type>()) L += two_pi;
-        if (L > mc::pi<promoted_type>()) L -= two_pi;
+        if (L < -pi) L += two_pi;
+        if (L > pi) L -= two_pi;
 
         if (math::equals(lat1, lat2) && math::equals(lon1, lon2))
         {
-            return return_type(0);
+            return calculation_type(0);
         }
 
         // TODO: give ellipsoid a template-parameter
-        promoted_type const ellipsoid_f = m_ellipsoid.f();
-        promoted_type const ellipsoid_a = m_ellipsoid.a();
-        promoted_type const ellipsoid_b = m_ellipsoid.b();
+        calculation_type const ellipsoid_f = m_ellipsoid.f();
+        calculation_type const ellipsoid_a = m_ellipsoid.a();
+        calculation_type const ellipsoid_b = m_ellipsoid.b();
 
         // U: reduced latitude, defined by tan U = (1-f) tan phi
-        promoted_type const c1 = 1;
-        promoted_type const one_min_f = c1 - ellipsoid_f;
+        calculation_type const c1 = 1;
+        calculation_type const one_min_f = c1 - ellipsoid_f;
 
-        promoted_type const U1 = atan(one_min_f * tan(lat1)); // above (1)
-        promoted_type const U2 = atan(one_min_f * tan(lat2)); // above (1)
+        calculation_type const U1 = atan(one_min_f * tan(lat1)); // above (1)
+        calculation_type const U2 = atan(one_min_f * tan(lat2)); // above (1)
 
-        promoted_type const cos_U1 = cos(U1);
-        promoted_type const cos_U2 = cos(U2);
-        promoted_type const sin_U1 = sin(U1);
-        promoted_type const sin_U2 = sin(U2);
+        calculation_type const cos_U1 = cos(U1);
+        calculation_type const cos_U2 = cos(U2);
+        calculation_type const sin_U1 = sin(U1);
+        calculation_type const sin_U2 = sin(U2);
 
         // alpha: azimuth of the geodesic at the equator
-        promoted_type cos2_alpha;
-        promoted_type sin_alpha;
+        calculation_type cos2_alpha;
+        calculation_type sin_alpha;
 
         // sigma: angular distance p1,p2 on the sphere
         // sigma1: angular distance on the sphere from the equator to p1
         // sigma_m: angular distance on the sphere from the equator to the midpoint of the line
-        promoted_type sigma;
-        promoted_type sin_sigma;
-        promoted_type cos2_sigma_m;
-
-        promoted_type previous_lambda;
-
-        promoted_type const c3 = 3;
-        promoted_type const c4 = 4;
-        promoted_type const c6 = 6;
-        promoted_type const c16 = 16;
+        calculation_type sigma;
+        calculation_type sin_sigma;
+        calculation_type cos2_sigma_m;
+
+        calculation_type previous_lambda;
+
+        calculation_type const c3 = 3;
+        calculation_type const c4 = 4;
+        calculation_type const c6 = 6;
+        calculation_type const c16 = 16;
 
-        promoted_type const c_e_12 = 1e-12;
+        calculation_type const c_e_12 = 1e-12;
 
         do
         {
             previous_lambda = lambda; // (13)
-            promoted_type sin_lambda = sin(lambda);
-            promoted_type cos_lambda = cos(lambda);
+            calculation_type sin_lambda = sin(lambda);
+            calculation_type cos_lambda = cos(lambda);
             sin_sigma = sqrt(math::sqr(cos_U2 * sin_lambda) + math::sqr(cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_lambda)); // (14)
-            promoted_type cos_sigma = sin_U1 * sin_U2 + cos_U1 * cos_U2 * cos_lambda; // (15)
+            calculation_type cos_sigma = sin_U1 * sin_U2 + cos_U1 * cos_U2 * cos_lambda; // (15)
             sin_alpha = cos_U1 * cos_U2 * sin_lambda / sin_sigma; // (17)
             cos2_alpha = c1 - math::sqr(sin_alpha);
-            cos2_sigma_m = cos2_alpha == 0 ? 0 : cos_sigma - c2 * sin_U1 * sin_U2 / cos2_alpha; // (18)
+            cos2_sigma_m = math::equals(cos2_alpha, 0) ? 0 : cos_sigma - c2 * sin_U1 * sin_U2 / cos2_alpha; // (18)
 
-
-            promoted_type C = ellipsoid_f/c16 * cos2_alpha * (c4 + ellipsoid_f * (c4 - c3 * cos2_alpha)); // (10)
+            calculation_type C = ellipsoid_f/c16 * cos2_alpha * (c4 + ellipsoid_f * (c4 - c3 * cos2_alpha)); // (10)
             sigma = atan2(sin_sigma, cos_sigma); // (16)
             lambda = L + (c1 - C) * ellipsoid_f * sin_alpha *
                 (sigma + C * sin_sigma * ( cos2_sigma_m + C * cos_sigma * (-c1 + c2 * math::sqr(cos2_sigma_m)))); // (11)
 
         } while (geometry::math::abs(previous_lambda - lambda) > c_e_12
-                && geometry::math::abs(lambda) < mc::pi<promoted_type>());
+                && geometry::math::abs(lambda) < pi);
 
-        promoted_type sqr_u = cos2_alpha * (math::sqr(ellipsoid_a) - math::sqr(ellipsoid_b)) / math::sqr(ellipsoid_b); // above (1)
+        calculation_type sqr_u = cos2_alpha * (math::sqr(ellipsoid_a) - math::sqr(ellipsoid_b)) / math::sqr(ellipsoid_b); // above (1)
 
         // Oops getting hard here
         // (again, problem is that ttmath cannot divide by doubles, which is OK)
-        promoted_type const c47 = 47;
-        promoted_type const c74 = 74;
-        promoted_type const c128 = 128;
-        promoted_type const c256 = 256;
-        promoted_type const c175 = 175;
-        promoted_type const c320 = 320;
-        promoted_type const c768 = 768;
-        promoted_type const c1024 = 1024;
-        promoted_type const c4096 = 4096;
-        promoted_type const c16384 = 16384;
-
-        promoted_type A = c1 + sqr_u/c16384 * (c4096 + sqr_u * (-c768 + sqr_u * (c320 - c175 * sqr_u))); // (3)
-        promoted_type B = sqr_u/c1024 * (c256 + sqr_u * ( -c128 + sqr_u * (c74 - c47 * sqr_u))); // (4)
-        promoted_type delta_sigma = B * sin_sigma * ( cos2_sigma_m + (B/c4) * (cos(sigma)* (-c1 + c2 * cos2_sigma_m)
+        calculation_type const c47 = 47;
+        calculation_type const c74 = 74;
+        calculation_type const c128 = 128;
+        calculation_type const c256 = 256;
+        calculation_type const c175 = 175;
+        calculation_type const c320 = 320;
+        calculation_type const c768 = 768;
+        calculation_type const c1024 = 1024;
+        calculation_type const c4096 = 4096;
+        calculation_type const c16384 = 16384;
+
+        calculation_type A = c1 + sqr_u/c16384 * (c4096 + sqr_u * (-c768 + sqr_u * (c320 - c175 * sqr_u))); // (3)
+        calculation_type B = sqr_u/c1024 * (c256 + sqr_u * ( -c128 + sqr_u * (c74 - c47 * sqr_u))); // (4)
+        calculation_type delta_sigma = B * sin_sigma * ( cos2_sigma_m + (B/c4) * (cos(sigma)* (-c1 + c2 * cos2_sigma_m)
                 - (B/c6) * cos2_sigma_m * (-c3 + c4 * math::sqr(sin_sigma)) * (-c3 + c4 * cos2_sigma_m))); // (6)
 
-        promoted_type dist = ellipsoid_b * A * (sigma - delta_sigma); // (19)
-
-        return return_type(dist);
+        return ellipsoid_b * A * (sigma - delta_sigma); // (19)
     }
 };
 
@@ -207,7 +202,7 @@
 template <typename Point1, typename Point2>
 struct return_type<strategy::distance::vincenty<Point1, Point2> >
 {
-    typedef typename strategy::distance::vincenty<Point1, Point2>::return_type type;
+    typedef typename strategy::distance::vincenty<Point1, Point2>::calculation_type type;
 };
 
 
@@ -247,7 +242,7 @@
 struct result_from_distance<vincenty<Point1, Point2> >
 {
     template <typename T>
-    static inline typename vincenty<Point1, Point2>::return_type apply(vincenty<Point1, Point2> const& , T const& value)
+    static inline typename return_type<vincenty<Point1, Point2> >::type apply(vincenty<Point1, Point2> const& , T const& value)
     {
         return value;
     }

Modified: sandbox/geometry/boost/geometry/util/math.hpp
==============================================================================
--- sandbox/geometry/boost/geometry/util/math.hpp	(original)
+++ sandbox/geometry/boost/geometry/util/math.hpp	2010-07-11 05:40:32 EDT (Sun, 11 Jul 2010)
@@ -19,9 +19,11 @@
 namespace boost { namespace geometry
 {
 
+namespace math
+{
 
 #ifndef DOXYGEN_NO_DETAIL
-namespace detail { namespace math
+namespace detail
 {
 
 
@@ -46,11 +48,26 @@
 };
 
 
-}} // namespace detail::math
+/*!
+\brief Short construct to enable partial specialization for PI, currently not possible in Math.
+*/
+template <typename T>
+struct define_pi
+{
+    static inline T apply()
+    {
+        // Default calls Boost.Math
+        return boost::math::constants::pi<T>();
+    }
+};
+
+
+} // namespace detail
 #endif
 
-namespace math
-{
+
+template <typename T>
+inline T pi() { return detail::define_pi<T>::apply(); }
 
 
 // Maybe replace this by boost equals or boost ublas numeric equals or so
@@ -72,7 +89,7 @@
 inline bool equals(T1 const& a, T2 const& b)
 {
     typedef typename select_most_precise<T1, T2>::type select_type;
-    return detail::math::equals
+    return detail::equals
         <
             select_type,
             boost::is_floating_point<select_type>::type::value
@@ -80,7 +97,7 @@
 }
 
 
-double const d2r = boost::math::constants::pi<double>() / 180.0;
+double const d2r = geometry::math::pi<double>() / 180.0;
 double const r2d = 1.0 / d2r;
 
 /*!
@@ -121,6 +138,7 @@
     return abs(t);
 }
 
+
 } // namespace math