|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r76755 - in trunk/boost/geometry: arithmetic policies/relate strategies/cartesian
From: barend.gehrels_at_[hidden]
Date: 2012-01-28 13:29:48
Author: barendgehrels
Date: 2012-01-28 13:29:47 EST (Sat, 28 Jan 2012)
New Revision: 76755
URL: http://svn.boost.org/trac/boost/changeset/76755
Log:
Introduced cross-product for area,centroid,side,intersection(determinant,direction,relation)
Text files modified:
trunk/boost/geometry/arithmetic/cross_product.hpp | 89 +++++++++++++++++++++++++++++++++++----
trunk/boost/geometry/policies/relate/direction.hpp | 30 +++++++++----
trunk/boost/geometry/policies/relate/intersection_points.hpp | 20 ++++----
trunk/boost/geometry/strategies/cartesian/area_surveyor.hpp | 11 +++-
trunk/boost/geometry/strategies/cartesian/cart_intersect.hpp | 13 ++++-
trunk/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp | 3
trunk/boost/geometry/strategies/cartesian/side_by_triangle.hpp | 11 +++-
7 files changed, 137 insertions(+), 40 deletions(-)
Modified: trunk/boost/geometry/arithmetic/cross_product.hpp
==============================================================================
--- trunk/boost/geometry/arithmetic/cross_product.hpp (original)
+++ trunk/boost/geometry/arithmetic/cross_product.hpp 2012-01-28 13:29:47 EST (Sat, 28 Jan 2012)
@@ -1,6 +1,8 @@
// Boost.Geometry (aka GGL, Generic Geometry Library)
// Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
+// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
@@ -12,8 +14,6 @@
#include <cstddef>
-#include <boost/concept/requires.hpp>
-
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/geometries/concepts/point_concept.hpp>
#include <boost/geometry/util/select_coordinate_type.hpp>
@@ -71,6 +71,30 @@
}
};
+
+template <typename ReturnType, typename U, typename V>
+class cross_product2
+{
+ template <typename T>
+ static inline ReturnType rt(T const& v)
+ {
+ return boost::numeric_cast<ReturnType>(v);
+ }
+
+public :
+
+ // Most common dimension, as also defined by Wolfram:
+ // http://mathworld.wolfram.com/CrossProduct.html
+ // "In two dimensions, the analog of the cross product for u=(u_x,u_y) and v=(v_x,v_y) is
+ // uxv = det(uv)
+ // = u_x v_y - u_y v_x"
+ static inline ReturnType apply(U const& ux, U const& uy
+ , V const& vx, V const& vy)
+ {
+ return rt(ux) * rt(vy) - rt(uy) * rt(vx);
+ }
+};
+
} // namespace detail
#endif // DOXYGEN_NO_DETAIL
@@ -83,13 +107,13 @@
// -- selection of lower dimension
/*!
- \brief Computes the cross product of two vector.
- \details Both vectors shall be of the same type.
- This type also determines type of result vector.
- \ingroup arithmetic
- \param p1 first vector
- \param p2 second vector
- \return the cross product vector
+\brief Computes the cross product of two vectors.
+\details Both vectors shall be of the same type.
+ This type also determines type of result vector.
+\ingroup arithmetic
+\param p1 first vector
+\param p2 second vector
+\return the cross product vector
*/
template <typename P1, typename P2>
inline P1 cross_product(P1 const& p1, P2 const& p2)
@@ -104,6 +128,53 @@
>::apply(p1, p2);
}
+
+/*!
+\brief Computes the cross product of two vectors, version for four values
+\details Because we often have the four coordinate values (often differences)
+ available, it is convenient to have a version which works directly on these,
+ without having to make a (temporary) Point or Vector
+\ingroup arithmetic
+\return the cross product value
+ */
+template <typename ReturnType, typename U, typename V>
+inline ReturnType cross_product2(U const& ux, U const& uy
+ , V const& vx, V const& vy)
+{
+ return detail::cross_product2
+ <
+ ReturnType, U, V
+ >::apply(ux, uy, vx, vy);
+}
+
+// Synonym, because yes, sometimes the algorithm calls it "determinant"
+template <typename ReturnType, typename U, typename V>
+inline ReturnType determinant(U const& ux, U const& uy
+ , V const& vx, V const& vy)
+{
+ return detail::cross_product2
+ <
+ ReturnType, U, V
+ >::apply(ux, uy, vx, vy);
+}
+
+
+// TEMPORARY, to be harmonized with cross_product
+template <typename ReturnType, typename U, typename V>
+inline ReturnType cross_product2(U const& u, V const& v)
+{
+ BOOST_CONCEPT_ASSERT( (concept::ConstPoint<U>) );
+ BOOST_CONCEPT_ASSERT( (concept::ConstPoint<V>) );
+
+ return detail::cross_product2
+ <
+ ReturnType,
+ typename geometry::coordinate_type<U>::type,
+ typename geometry::coordinate_type<V>::type
+ >::apply(get<0>(u), get<1>(u), get<0>(v), get<1>(v));
+}
+
+
}} // namespace boost::geometry
#endif // BOOST_GEOMETRY_ARITHMETIC_CROSS_PRODUCT_HPP
Modified: trunk/boost/geometry/policies/relate/direction.hpp
==============================================================================
--- trunk/boost/geometry/policies/relate/direction.hpp (original)
+++ trunk/boost/geometry/policies/relate/direction.hpp 2012-01-28 13:29:47 EST (Sat, 28 Jan 2012)
@@ -15,6 +15,7 @@
#include <boost/concept_check.hpp>
+#include <boost/geometry/arithmetic/cross_product.hpp>
#include <boost/geometry/strategies/side_info.hpp>
#include <boost/geometry/util/math.hpp>
@@ -249,6 +250,21 @@
private :
+ static inline bool is_left
+ (
+ coordinate_type const& ux,
+ coordinate_type const& uy,
+ coordinate_type const& vx,
+ coordinate_type const& vy
+ )
+ {
+ // This is a "side calculation" as in the strategies, but here terms are precalculated
+ // We might merge this with side, offering a pre-calculated term (in fact already done using cross-product)
+ // Waiting for implementing spherical...
+
+ rtype const zero = rtype();
+ return geometry::cross_product2<rtype>(ux, uy, vx, vy) > zero;
+ }
template <std::size_t I>
static inline return_type calculate_side(side_info const& sides,
@@ -259,11 +275,7 @@
coordinate_type dpx = get<I, 0>(s2) - get<0, 0>(s1);
coordinate_type dpy = get<I, 1>(s2) - get<0, 1>(s1);
- // This is a "side calculation" as in the strategies, but here two terms are precalculated
- // We might merge this with side, offering a pre-calculated term
- // Waiting for implementing spherical...
-
- return dx1 * dpy - dy1 * dpx > 0
+ return is_left(dx1, dy1, dpx, dpy)
? return_type(sides, how, how_a, how_b, -1, 1)
: return_type(sides, how, how_a, how_b, 1, -1);
}
@@ -277,7 +289,7 @@
coordinate_type dpx = get<I, 0>(s2) - get<0, 0>(s1);
coordinate_type dpy = get<I, 1>(s2) - get<0, 1>(s1);
- return dx1 * dpy - dy1 * dpx > 0
+ return is_left(dx1, dy1, dpx, dpy)
? return_type(sides, how, how_a, how_b, 1, 1)
: return_type(sides, how, how_a, how_b, -1, -1);
}
@@ -293,7 +305,7 @@
coordinate_type dpx = get<1, 0>(s2) - get<0, 0>(s1);
coordinate_type dpy = get<1, 1>(s2) - get<0, 1>(s1);
- int dir = dx1 * dpy - dy1 * dpx > 0 ? 1 : -1;
+ int dir = is_left(dx1, dy1, dpx, dpy) ? 1 : -1;
// From other perspective, then reverse
bool const is_a = which == 'A';
@@ -321,7 +333,7 @@
// Ending at the middle, one ARRIVES, the other one is NEUTRAL
// (because it both "arrives" and "departs" there
- return dx * dpy - dy * dpx > 0
+ return is_left(dx, dy, dpx, dpy)
? return_type(sides, 'm', 1, 0, 1, 1)
: return_type(sides, 'm', 1, 0, -1, -1);
}
@@ -334,7 +346,7 @@
coordinate_type dpx = get<1, 0>(s1) - get<0, 0>(s2);
coordinate_type dpy = get<1, 1>(s1) - get<0, 1>(s2);
- return dx * dpy - dy * dpx > 0
+ return is_left(dx, dy, dpx, dpy)
? return_type(sides, 'm', 0, 1, 1, 1)
: return_type(sides, 'm', 0, 1, -1, -1);
}
Modified: trunk/boost/geometry/policies/relate/intersection_points.hpp
==============================================================================
--- trunk/boost/geometry/policies/relate/intersection_points.hpp (original)
+++ trunk/boost/geometry/policies/relate/intersection_points.hpp 2012-01-28 13:29:47 EST (Sat, 28 Jan 2012)
@@ -16,6 +16,7 @@
#include <boost/concept_check.hpp>
#include <boost/numeric/conversion/cast.hpp>
+#include <boost/geometry/arithmetic/cross_product.hpp>
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/strategies/side_info.hpp>
#include <boost/geometry/util/select_calculation_type.hpp>
@@ -40,9 +41,6 @@
S1, S2, CalculationType
>::type coordinate_type;
- // Get the same type, but at least a double
- typedef typename select_most_precise<coordinate_type, double>::type rtype;
-
static inline return_type segments_intersect(side_info const&,
coordinate_type const& dx1, coordinate_type const& dy1,
coordinate_type const& dx2, coordinate_type const& dy2,
@@ -54,7 +52,7 @@
typename return_type::point_type
>::type coordinate_type;
- // Get the same type, but at least a double (also used for divisions
+ // Get the same type, but at least a double (also used for divisions)
typedef typename select_most_precise
<
coordinate_type, double
@@ -66,19 +64,21 @@
// Calculate other determinants - Cramers rule
promoted_type const wx = get<0, 0>(s1) - get<0, 0>(s2);
promoted_type const wy = get<0, 1>(s1) - get<0, 1>(s2);
- promoted_type const d = (dy2 * dx1) - (dx2 * dy1);
- promoted_type const da = (promoted_type(dx2) * wy) - (promoted_type(dy2) * wx);
+ promoted_type const d = geometry::determinant<promoted_type>(dx1, dy1, dx2, dy2);
+ promoted_type const da = geometry::determinant<promoted_type>(dx2, dy2, wx, wy);
// r: ratio 0-1 where intersection divides A/B
promoted_type r = da / d;
+ promoted_type const zero = promoted_type();
+ promoted_type const one = 1;
// Handle robustness issues
- if (r < 0)
+ if (r < zero)
{
- r = 0;
+ r = zero;
}
- else if (r > 1)
+ else if (r > one)
{
- r = 1;
+ r = one;
}
result.count = 1;
Modified: trunk/boost/geometry/strategies/cartesian/area_surveyor.hpp
==============================================================================
--- trunk/boost/geometry/strategies/cartesian/area_surveyor.hpp (original)
+++ trunk/boost/geometry/strategies/cartesian/area_surveyor.hpp 2012-01-28 13:29:47 EST (Sat, 28 Jan 2012)
@@ -16,9 +16,11 @@
#include <boost/mpl/if.hpp>
-#include <boost/type_traits.hpp>
-#include <boost/geometry/geometries/point_xy.hpp>
+#include <boost/geometry/arithmetic/cross_product.hpp>
+#include <boost/geometry/core/coordinate_type.hpp>
+#include <boost/geometry/core/coordinate_dimension.hpp>
+#include <boost/geometry/util/select_most_precise.hpp>
namespace boost { namespace geometry
@@ -82,7 +84,8 @@
inline return_type area() const
{
return_type result = sum;
- result /= 2;
+ return_type const two = 2;
+ result /= two;
return result;
}
};
@@ -96,7 +99,7 @@
summation& state)
{
// SUM += x2 * y1 - x1 * y2;
- state.sum += get<0>(p2) * get<1>(p1) - get<0>(p1) * get<1>(p2);
+ state.sum += geometry::cross_product2<return_type>(p2, p1);
}
static inline return_type result(summation const& state)
Modified: trunk/boost/geometry/strategies/cartesian/cart_intersect.hpp
==============================================================================
--- trunk/boost/geometry/strategies/cartesian/cart_intersect.hpp (original)
+++ trunk/boost/geometry/strategies/cartesian/cart_intersect.hpp 2012-01-28 13:29:47 EST (Sat, 28 Jan 2012)
@@ -16,6 +16,7 @@
#include <boost/geometry/geometries/concepts/point_concept.hpp>
#include <boost/geometry/geometries/concepts/segment_concept.hpp>
+#include <boost/geometry/arithmetic/cross_product.hpp>
#include <boost/geometry/algorithms/detail/assign_values.hpp>
#include <boost/geometry/util/math.hpp>
@@ -199,13 +200,19 @@
coordinate_type, double
>::type promoted_type;
+ // Calculate the determinant/2D cross product
+ // (Note, because we only check on zero,
+ // the order a/b does not care)
+ promoted_type const d = geometry::determinant
+ <
+ promoted_type
+ >(dx_a, dy_a, dx_b, dy_b);
- promoted_type const d = (dy_b * dx_a) - (dx_b * dy_a);
// Determinant d should be nonzero.
// If it is zero, we have an robustness issue here,
// (and besides that we cannot divide by it)
- if(math::equals(d, zero) && ! collinear)
- //if(! collinear && sides.as_collinear())
+ promoted_type const pt_zero = promoted_type();
+ if(math::equals(d, pt_zero) && ! collinear)
{
#ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION
std::cout << "Determinant zero? Type : "
Modified: trunk/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp
==============================================================================
--- trunk/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp (original)
+++ trunk/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp 2012-01-28 13:29:47 EST (Sat, 28 Jan 2012)
@@ -19,6 +19,7 @@
#include <boost/numeric/conversion/cast.hpp>
#include <boost/type_traits.hpp>
+#include <boost/geometry/arithmetic/cross_product.hpp>
#include <boost/geometry/core/coordinate_type.hpp>
#include <boost/geometry/core/point_type.hpp>
#include <boost/geometry/strategies/centroid.hpp>
@@ -177,7 +178,7 @@
calculation_type const y1 = boost::numeric_cast<calculation_type>(get<1>(p1));
calculation_type const x2 = boost::numeric_cast<calculation_type>(get<0>(p2));
calculation_type const y2 = boost::numeric_cast<calculation_type>(get<1>(p2));
- calculation_type const ai = x1 * y2 - x2 * y1;
+ calculation_type const ai = geometry::cross_product2<calculation_type>(p1, p2);
state.count++;
state.sum_a2 += ai;
state.sum_x += ai * (x1 + x2);
Modified: trunk/boost/geometry/strategies/cartesian/side_by_triangle.hpp
==============================================================================
--- trunk/boost/geometry/strategies/cartesian/side_by_triangle.hpp (original)
+++ trunk/boost/geometry/strategies/cartesian/side_by_triangle.hpp 2012-01-28 13:29:47 EST (Sat, 28 Jan 2012)
@@ -17,10 +17,9 @@
#include <boost/mpl/if.hpp>
#include <boost/type_traits.hpp>
+#include <boost/geometry/arithmetic/cross_product.hpp>
#include <boost/geometry/core/access.hpp>
-
#include <boost/geometry/util/select_coordinate_type.hpp>
-
#include <boost/geometry/strategies/side.hpp>
@@ -66,7 +65,6 @@
CalculationType
>::type coordinate_type;
-//std::cout << "side: " << typeid(coordinate_type).name() << std::endl;
coordinate_type const x = get<0>(p);
coordinate_type const y = get<1>(p);
@@ -87,7 +85,12 @@
promoted_type const dpx = x - sx1;
promoted_type const dpy = y - sy1;
- promoted_type const s = dx * dpy - dy * dpx;
+ promoted_type const s
+ = geometry::cross_product2<promoted_type>
+ (
+ dx, dy,
+ dpx, dpy
+ );
promoted_type const zero = promoted_type();
return math::equals(s, zero) ? 0
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