|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r63842 - in sandbox/geometry: boost/geometry/extensions/gis/geographic/detail boost/geometry/extensions/gis/geographic/strategies libs/geometry/test/extensions/gis/latlong
From: barend.gehrels_at_[hidden]
Date: 2010-07-11 06:57:48
Author: barendgehrels
Date: 2010-07-11 06:57:45 EDT (Sun, 11 Jul 2010)
New Revision: 63842
URL: http://svn.boost.org/trac/boost/changeset/63842
Log:
Added template parameter to ellipsoid
Text files modified:
sandbox/geometry/boost/geometry/extensions/gis/geographic/detail/ellipsoid.hpp | 17 ++++++------
sandbox/geometry/boost/geometry/extensions/gis/geographic/strategies/andoyer.hpp | 53 ++++++++++++++++++---------------------
sandbox/geometry/boost/geometry/extensions/gis/geographic/strategies/vincenty.hpp | 21 ++++++---------
sandbox/geometry/libs/geometry/test/extensions/gis/latlong/andoyer.cpp | 13 +++++++++
sandbox/geometry/libs/geometry/test/extensions/gis/latlong/andoyer.vcproj | 8 +++---
sandbox/geometry/libs/geometry/test/extensions/gis/latlong/vincenty.cpp | 5 ++-
6 files changed, 61 insertions(+), 56 deletions(-)
Modified: sandbox/geometry/boost/geometry/extensions/gis/geographic/detail/ellipsoid.hpp
==============================================================================
--- sandbox/geometry/boost/geometry/extensions/gis/geographic/detail/ellipsoid.hpp (original)
+++ sandbox/geometry/boost/geometry/extensions/gis/geographic/detail/ellipsoid.hpp 2010-07-11 06:57:45 EDT (Sun, 11 Jul 2010)
@@ -23,31 +23,32 @@
- and http://en.wikipedia.org/wiki/World_Geodetic_System#A_new_World_Geodetic_System:_WGS84
\note
*/
+template <typename T>
class ellipsoid
{
public :
- ellipsoid(double a, double b)
+ ellipsoid(T const& a, T const& b)
: m_a(a)
, m_b(b)
, m_f((a - b) / a)
{}
ellipsoid()
- : m_a(6378137.0)
- , m_b(6356752.314245)
+ : m_a(T(6378137.0))
+ , m_b(T(6356752.314245))
, m_f((m_a - m_b) / m_a)
{}
// Unit sphere
- ellipsoid(double f)
+ ellipsoid(T const& f)
: m_a(1.0)
, m_f(f)
{}
- double a() const { return m_a; }
- double b() const { return m_b; }
- double f() const { return m_f; }
+ T a() const { return m_a; }
+ T b() const { return m_b; }
+ T f() const { return m_f; }
private :
- double m_a, m_b, m_f; // equatorial radius, polar radius, flattening
+ T m_a, m_b, m_f; // equatorial radius, polar radius, flattening
};
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 06:57:45 EDT (Sun, 11 Jul 2010)
@@ -69,7 +69,7 @@
: m_ellipsoid(f)
{}
- explicit inline andoyer(geometry::detail::ellipsoid const& e)
+ explicit inline andoyer(geometry::detail::ellipsoid<calculation_type> const& e)
: m_ellipsoid(e)
{}
@@ -80,24 +80,22 @@
get_as_radian<0>(point2), get_as_radian<1>(point2));
}
- inline geometry::detail::ellipsoid ellipsoid() const
+ inline geometry::detail::ellipsoid<calculation_type> ellipsoid() const
{
return m_ellipsoid;
}
private :
- typedef typename coordinate_type<Point1>::type T1;
- typedef typename coordinate_type<Point2>::type T2;
- geometry::detail::ellipsoid m_ellipsoid;
+ geometry::detail::ellipsoid<calculation_type> m_ellipsoid;
inline calculation_type calc(calculation_type const& lon1,
calculation_type const& lat1,
calculation_type const& lon2,
calculation_type const& lat2) const
{
- calculation_type G = (lat1 - lat2) / 2.0;
- calculation_type lambda = (lon1 - lon2) / 2.0;
+ calculation_type const G = (lat1 - lat2) / 2.0;
+ calculation_type const lambda = (lon1 - lon2) / 2.0;
if (geometry::math::equals(lambda, 0.0)
&& geometry::math::equals(G, 0.0))
@@ -105,35 +103,34 @@
return 0.0;
}
- calculation_type F = (lat1 + lat2) / 2.0;
+ calculation_type const F = (lat1 + lat2) / 2.0;
- calculation_type sinG2 = math::sqr(sin(G));
- calculation_type cosG2 = math::sqr(cos(G));
- calculation_type sinF2 = math::sqr(sin(F));
- calculation_type cosF2 = math::sqr(cos(F));
- calculation_type sinL2 = math::sqr(sin(lambda));
- calculation_type cosL2 = math::sqr(cos(lambda));
+ calculation_type const sinG2 = math::sqr(sin(G));
+ calculation_type const cosG2 = math::sqr(cos(G));
+ calculation_type const sinF2 = math::sqr(sin(F));
+ calculation_type const cosF2 = math::sqr(cos(F));
+ calculation_type const sinL2 = math::sqr(sin(lambda));
+ calculation_type const cosL2 = math::sqr(cos(lambda));
- calculation_type S = sinG2 * cosL2 + cosF2 * sinL2;
- calculation_type C = cosG2 * cosL2 + sinF2 * sinL2;
-
- if (geometry::math::equals(S, 0.0) || geometry::math::equals(C, 0.0))
- {
- return 0.0;
- }
+ calculation_type const S = sinG2 * cosL2 + cosF2 * sinL2;
+ calculation_type const C = cosG2 * cosL2 + sinF2 * sinL2;
+ calculation_type const c0 = 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 = c2 * omega * m_ellipsoid.a();
- calculation_type H1 = (c3 * r - c1) / (c2 * C);
- calculation_type H2 = (c3 * r + c1) / (c2 * S);
+ if (geometry::math::equals(S, c0) || geometry::math::equals(C, c0))
+ {
+ return c0;
+ }
- calculation_type f = m_ellipsoid.f();
+ calculation_type const omega = atan(sqrt(S / C));
+ calculation_type const r3 = c3 * sqrt(S * C) / omega; // not sure if this is r or greek nu
+ calculation_type const D = c2 * omega * m_ellipsoid.a();
+ calculation_type const H1 = (r3 - c1) / (c2 * C);
+ calculation_type const H2 = (r3 + c1) / (c2 * S);
+ calculation_type const f = m_ellipsoid.f();
return D * (c1 + f * H1 * sinF2 * cosG2 - f * H2 * cosF2 * sinG2);
}
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 06:57:45 EDT (Sun, 11 Jul 2010)
@@ -68,7 +68,7 @@
inline vincenty()
{}
- explicit inline vincenty(geometry::detail::ellipsoid const& e)
+ explicit inline vincenty(geometry::detail::ellipsoid<calculation_type> const& e)
: m_ellipsoid(e)
{}
@@ -78,14 +78,14 @@
get_as_radian<0>(p2), get_as_radian<1>(p2));
}
- inline geometry::detail::ellipsoid ellipsoid() const
+ inline geometry::detail::ellipsoid<calculation_type> ellipsoid() const
{
return m_ellipsoid;
}
private :
- geometry::detail::ellipsoid m_ellipsoid;
+ geometry::detail::ellipsoid<calculation_type> m_ellipsoid;
inline calculation_type calculate(calculation_type const& lon1,
calculation_type const& lat1,
@@ -108,14 +108,9 @@
return calculation_type(0);
}
- // TODO: give ellipsoid a template-parameter
- 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
calculation_type const c1 = 1;
- calculation_type const one_min_f = c1 - ellipsoid_f;
+ calculation_type const one_min_f = c1 - m_ellipsoid.f();
calculation_type const U1 = atan(one_min_f * tan(lat1)); // above (1)
calculation_type const U2 = atan(one_min_f * tan(lat2)); // above (1)
@@ -156,15 +151,15 @@
cos2_alpha = c1 - math::sqr(sin_alpha);
cos2_sigma_m = math::equals(cos2_alpha, 0) ? 0 : cos_sigma - c2 * sin_U1 * sin_U2 / cos2_alpha; // (18)
- calculation_type C = ellipsoid_f/c16 * cos2_alpha * (c4 + ellipsoid_f * (c4 - c3 * cos2_alpha)); // (10)
+ calculation_type C = m_ellipsoid.f()/c16 * cos2_alpha * (c4 + m_ellipsoid.f() * (c4 - c3 * cos2_alpha)); // (10)
sigma = atan2(sin_sigma, cos_sigma); // (16)
- lambda = L + (c1 - C) * ellipsoid_f * sin_alpha *
+ lambda = L + (c1 - C) * m_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) < pi);
- calculation_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(m_ellipsoid.a()) - math::sqr(m_ellipsoid.b())) / math::sqr(m_ellipsoid.b()); // above (1)
// Oops getting hard here
// (again, problem is that ttmath cannot divide by doubles, which is OK)
@@ -184,7 +179,7 @@
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)
- return ellipsoid_b * A * (sigma - delta_sigma); // (19)
+ return m_ellipsoid.b() * A * (sigma - delta_sigma); // (19)
}
};
Modified: sandbox/geometry/libs/geometry/test/extensions/gis/latlong/andoyer.cpp
==============================================================================
--- sandbox/geometry/libs/geometry/test/extensions/gis/latlong/andoyer.cpp (original)
+++ sandbox/geometry/libs/geometry/test/extensions/gis/latlong/andoyer.cpp 2010-07-11 06:57:45 EDT (Sun, 11 Jul 2010)
@@ -18,6 +18,10 @@
#include <boost/geometry/geometries/point.hpp>
#include <test_common/test_point.hpp>
+#ifdef HAVE_TTMATH
+# include <boost/geometry/extensions/contrib/ttmath_stub.hpp>
+#endif
+
namespace bg = boost::geometry;
@@ -29,6 +33,7 @@
BOOST_CONCEPT_ASSERT( (bg::concept::PointDistanceStrategy<andoyer_type>) );
andoyer_type andoyer;
+ typedef bg::strategy::distance::services::return_type<andoyer_type>::type return_type;
P1 p1, p2;
@@ -36,7 +41,7 @@
bg::assign(p1, lon1, lat1);
bg::assign(p2, lon2, lat2);
- BOOST_CHECK_CLOSE(double(andoyer.apply(p1, p2)), 1000.0 * expected_km, 0.001);
+ BOOST_CHECK_CLOSE(andoyer.apply(p1, p2), return_type(1000.0 * expected_km), 0.001);
}
template <typename P1, typename P2>
@@ -57,8 +62,14 @@
{
//test_all<float[2]>();
//test_all<double[2]>();
+ test_all<bg::point<int, 2, bg::cs::geographic<bg::degree> > >();
test_all<bg::point<float, 2, bg::cs::geographic<bg::degree> > >();
test_all<bg::point<double, 2, bg::cs::geographic<bg::degree> > >();
+#if defined(HAVE_TTMATH)
+ test_all<bg::point<ttmath::Big<1,4>, 2, bg::cs::geographic<bg::degree> > >();
+ test_all<bg::point<ttmath_big, 2, bg::cs::geographic<bg::degree> > >();
+#endif
+
return 0;
}
Modified: sandbox/geometry/libs/geometry/test/extensions/gis/latlong/andoyer.vcproj
==============================================================================
--- sandbox/geometry/libs/geometry/test/extensions/gis/latlong/andoyer.vcproj (original)
+++ sandbox/geometry/libs/geometry/test/extensions/gis/latlong/andoyer.vcproj 2010-07-11 06:57:45 EDT (Sun, 11 Jul 2010)
@@ -20,7 +20,7 @@
OutputDirectory="$(SolutionDir)$(ConfigurationName)"
IntermediateDirectory="$(ConfigurationName)\andoyer"
ConfigurationType="1"
- InheritedPropertySheets="..\..\..\boost.vsprops"
+ InheritedPropertySheets="..\..\..\boost.vsprops;..\..\..\ttmath.vsprops"
CharacterSet="1"
>
<Tool
@@ -41,7 +41,7 @@
<Tool
Name="VCCLCompilerTool"
Optimization="0"
- AdditionalIncludeDirectories="../../../../../..;../../.."
+ AdditionalIncludeDirectories="../../../../../..;../../..;../../../../../../boost/geometry/extensions/contrib"
PreprocessorDefinitions="WIN32;_DEBUG;_CONSOLE"
ExceptionHandling="2"
RuntimeLibrary="1"
@@ -93,7 +93,7 @@
OutputDirectory="$(SolutionDir)$(ConfigurationName)"
IntermediateDirectory="$(ConfigurationName)\andoyer"
ConfigurationType="1"
- InheritedPropertySheets="..\..\..\boost.vsprops"
+ InheritedPropertySheets="..\..\..\boost.vsprops;..\..\..\ttmath.vsprops"
CharacterSet="1"
WholeProgramOptimization="1"
>
@@ -114,7 +114,7 @@
/>
<Tool
Name="VCCLCompilerTool"
- AdditionalIncludeDirectories="../../../../../..;../../.."
+ AdditionalIncludeDirectories="../../../../../..;../../..;../../../../../../boost/geometry/extensions/contrib"
PreprocessorDefinitions="WIN32;NDEBUG;_CONSOLE"
ExceptionHandling="2"
UsePrecompiledHeader="0"
Modified: sandbox/geometry/libs/geometry/test/extensions/gis/latlong/vincenty.cpp
==============================================================================
--- sandbox/geometry/libs/geometry/test/extensions/gis/latlong/vincenty.cpp (original)
+++ sandbox/geometry/libs/geometry/test/extensions/gis/latlong/vincenty.cpp 2010-07-11 06:57:45 EDT (Sun, 11 Jul 2010)
@@ -24,6 +24,7 @@
namespace bg = boost::geometry;
+
template <typename P1, typename P2>
void test_vincenty(double lon1, double lat1, double lon2, double lat2, double expected_km)
{
@@ -46,7 +47,7 @@
template <typename P1, typename P2>
void test_all()
{
- test_vincenty<P1, P2>(0, 90, 1, 80, 1116.828862); // polar
+ test_vincenty<P1, P2>(0, 89, 1, 80, 1005.1535769); // sub-polar
test_vincenty<P1, P2>(4, 52, 4, 52, 0.0); // no point difference
test_vincenty<P1, P2>(4, 52, 3, 40, 1336.039890); // normal case
}
@@ -63,7 +64,7 @@
//test_all<float[2]>();
//test_all<double[2]>();
test_all<bg::point<int, 2, bg::cs::geographic<bg::degree> > >();
- //test_all<point<float, 2, geographic<degree> > >();
+ test_all<bg::point<float, 2, bg::cs::geographic<bg::degree> > >();
test_all<bg::point<double, 2, bg::cs::geographic<bg::degree> > >();
#if defined(HAVE_TTMATH)
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk