Subject: [Boost-commit] svn:boost r72286 - in trunk/boost/geometry/strategies: cartesian spherical
From: barend.gehrels_at_[hidden]
Date: 2011-05-30 11:07:13
Author: barendgehrels
Date: 2011-05-30 11:07:12 EDT (Mon, 30 May 2011)
New Revision: 72286
URL: http://svn.boost.org/trac/boost/changeset/72286
Log:
Huiller: changed to calculation types (now supporting ttmath)
distance_projected_point.hpp: minor changes
Text files modified:
trunk/boost/geometry/strategies/cartesian/distance_projected_point.hpp | 25 ++++++++++-----------
trunk/boost/geometry/strategies/spherical/area_huiller.hpp | 47 +++++++++++++++++++++++----------------
2 files changed, 40 insertions(+), 32 deletions(-)
Modified: trunk/boost/geometry/strategies/cartesian/distance_projected_point.hpp
==============================================================================
--- trunk/boost/geometry/strategies/cartesian/distance_projected_point.hpp (original)
+++ trunk/boost/geometry/strategies/cartesian/distance_projected_point.hpp 2011-05-30 11:07:12 EDT (Mon, 30 May 2011)
@@ -116,13 +116,14 @@
{
assert_dimension_equal<Point, PointOfSegment>();
- /* Algorithm
- POINT v(x2 - x1, y2 - y1);
- POINT w(px - x1, py - y1);
- c1 = w . v
- c2 = v . v
- b = c1 / c2
- RETURN POINT(x1 + b * vx, y1 + b * vy);
+ /*
+ Algorithm [p1: (x1,y1), p2: (x2,y2), p: (px,py)]
+ VECTOR v(x2 - x1, y2 - y1)
+ VECTOR w(px - x1, py - y1)
+ c1 = w . v
+ c2 = v . v
+ b = c1 / c2
+ RETURN POINT(x1 + b * vx, y1 + b * vy)
*/
// v is multiplied below with a (possibly) FP-value, so should be in FP
@@ -137,21 +138,20 @@
Strategy strategy;
boost::ignore_unused_variable_warning(strategy);
- calculation_type zero = calculation_type();
- fp_type c1 = dot_product(w, v);
+ calculation_type const zero = calculation_type();
+ fp_type const c1 = dot_product(w, v);
if (c1 <= zero)
{
return strategy.apply(p, p1);
}
- fp_type c2 = dot_product(v, v);
+ fp_type const c2 = dot_product(v, v);
if (c2 <= c1)
{
return strategy.apply(p, p2);
}
// See above, c1 > 0 AND c2 > c1 so: c2 != 0
- fp_type b = fp_type(c1) / fp_type(c2);
-
+ fp_type const b = c1 / c2;
fp_strategy_type fp_strategy
= strategy::distance::services::get_similar
@@ -168,7 +168,6 @@
return fp_strategy.apply(p, projected);
}
-
};
#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
Modified: trunk/boost/geometry/strategies/spherical/area_huiller.hpp
==============================================================================
--- trunk/boost/geometry/strategies/spherical/area_huiller.hpp (original)
+++ trunk/boost/geometry/strategies/spherical/area_huiller.hpp 2011-05-30 11:07:12 EDT (Mon, 30 May 2011)
@@ -10,7 +10,6 @@
#define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_HUILLER_HPP
-#include <boost/math/constants/constants.hpp>
#include <boost/geometry/strategies/spherical/distance_haversine.hpp>
@@ -66,10 +65,21 @@
>
class huiller
{
+typedef typename boost::mpl::if_c
+ <
+ boost::is_void<CalculationType>::type::value,
+ typename select_most_precise
+ <
+ typename coordinate_type<PointOfSegment>::type,
+ double
+ >::type,
+ CalculationType
+ >::type calculation_type;
+
protected :
struct excess_sum
{
- double sum;
+ calculation_type sum;
// Distances are calculated on unit sphere here
strategy::distance::haversine<PointOfSegment, PointOfSegment>
@@ -80,18 +90,18 @@
: sum(0)
, distance_over_unit_sphere(1)
{}
- inline double area(double radius) const
+ inline calculation_type area(calculation_type radius) const
{
return - sum * radius * radius;
}
};
public :
- typedef double return_type;
+ typedef calculation_type return_type;
typedef PointOfSegment segment_point_type;
typedef excess_sum state_type;
- inline huiller(double radius = 1.0)
+ inline huiller(calculation_type radius = 1.0)
: m_radius(radius)
{}
@@ -101,26 +111,25 @@
{
if (! geometry::math::equals(get<0>(p1), get<0>(p2)))
{
- namespace mc = boost::math::constants;
-
- double const two_pi = 2.0 * mc::pi<double>();
- double const half = 0.5;
- double const two = 2.0;
- double const four = 4.0;
+ calculation_type const half = 0.5;
+ calculation_type const two = 2.0;
+ calculation_type const four = 4.0;
+ calculation_type const two_pi = two * geometry::math::pi<calculation_type>();
+ calculation_type const half_pi = half * geometry::math::pi<calculation_type>();
// Distance p1 p2
- double a = state.distance_over_unit_sphere.apply(p1, p2);
+ calculation_type a = state.distance_over_unit_sphere.apply(p1, p2);
// Sides on unit sphere to south pole
- double b = half * mc::pi<double>() - geometry::get_as_radian<1>(p2);
- double c = half * mc::pi<double>() - geometry::get_as_radian<1>(p1);
+ calculation_type b = half_pi - geometry::get_as_radian<1>(p2);
+ calculation_type c = half_pi - geometry::get_as_radian<1>(p1);
// Semi parameter
- double s = half * (a + b + c);
+ calculation_type s = half * (a + b + c);
// E: spherical excess, using l'Huiller's formula
// [tg(e / 4)]2 = tg[s / 2] tg[(s-a) / 2] tg[(s-b) / 2] tg[(s-c) / 2]
- double E = four * atan(sqrt(geometry::math::abs(tan(s / two)
+ calculation_type E = four * atan(sqrt(geometry::math::abs(tan(s / two)
* tan((s - a) / two)
* tan((s - b) / two)
* tan((s - c) / two))));
@@ -132,11 +141,11 @@
// we have to take the dateline into account.
// TODO: check this / enhance this, should be more robust. See also the "grow" for ll
// TODO: use minmax or "smaller"/"compare" strategy for this
- double lon1 = geometry::get_as_radian<0>(p1) < 0
+ calculation_type lon1 = geometry::get_as_radian<0>(p1) < 0
? geometry::get_as_radian<0>(p1) + two_pi
: geometry::get_as_radian<0>(p1);
- double lon2 = geometry::get_as_radian<0>(p2) < 0
+ calculation_type lon2 = geometry::get_as_radian<0>(p2) < 0
? geometry::get_as_radian<0>(p2) + two_pi
: geometry::get_as_radian<0>(p2);
@@ -156,7 +165,7 @@
private :
/// Radius of the sphere
- double m_radius;
+ calculation_type m_radius;
};
#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS