From 4594a65da420533205df3dcc3ae0db05ce91c008 Mon Sep 17 00:00:00 2001 From: Barend Gehrels Date: Sat, 28 Jan 2012 18:29:47 +0000 Subject: [PATCH] Introduced cross-product for area,centroid,side,intersection(determinant,direction,relation) [SVN r76755] --- .../geometry/arithmetic/cross_product.hpp | 89 +++++++++++++++++-- .../geometry/policies/relate/direction.hpp | 30 +++++-- .../policies/relate/intersection_points.hpp | 20 ++--- .../strategies/cartesian/area_surveyor.hpp | 11 ++- .../strategies/cartesian/cart_intersect.hpp | 13 ++- .../cartesian/centroid_bashein_detmer.hpp | 3 +- .../strategies/cartesian/side_by_triangle.hpp | 11 ++- 7 files changed, 137 insertions(+), 40 deletions(-) diff --git a/include/boost/geometry/arithmetic/cross_product.hpp b/include/boost/geometry/arithmetic/cross_product.hpp index e7ec1a59b..da90ce40d 100644 --- a/include/boost/geometry/arithmetic/cross_product.hpp +++ b/include/boost/geometry/arithmetic/cross_product.hpp @@ -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 -#include - #include #include #include @@ -71,6 +71,30 @@ struct cross_product } }; + +template +class cross_product2 +{ + template + static inline ReturnType rt(T const& v) + { + return boost::numeric_cast(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 @@ struct cross_product // -- 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 inline P1 cross_product(P1 const& p1, P2 const& p2) @@ -104,6 +128,53 @@ inline P1 cross_product(P1 const& p1, P2 const& p2) >::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 +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 +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 +inline ReturnType cross_product2(U const& u, V const& v) +{ + BOOST_CONCEPT_ASSERT( (concept::ConstPoint) ); + BOOST_CONCEPT_ASSERT( (concept::ConstPoint) ); + + return detail::cross_product2 + < + ReturnType, + typename geometry::coordinate_type::type, + typename geometry::coordinate_type::type + >::apply(get<0>(u), get<1>(u), get<0>(v), get<1>(v)); +} + + }} // namespace boost::geometry #endif // BOOST_GEOMETRY_ARITHMETIC_CROSS_PRODUCT_HPP diff --git a/include/boost/geometry/policies/relate/direction.hpp b/include/boost/geometry/policies/relate/direction.hpp index ba190fa04..34c8b668b 100644 --- a/include/boost/geometry/policies/relate/direction.hpp +++ b/include/boost/geometry/policies/relate/direction.hpp @@ -15,6 +15,7 @@ #include +#include #include #include @@ -249,6 +250,21 @@ struct segments_direction 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(ux, uy, vx, vy) > zero; + } template static inline return_type calculate_side(side_info const& sides, @@ -259,11 +275,7 @@ private : coordinate_type dpx = get(s2) - get<0, 0>(s1); coordinate_type dpy = get(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 @@ private : coordinate_type dpx = get(s2) - get<0, 0>(s1); coordinate_type dpy = get(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 @@ private : 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 @@ private : // 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 @@ private : 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); } diff --git a/include/boost/geometry/policies/relate/intersection_points.hpp b/include/boost/geometry/policies/relate/intersection_points.hpp index 5d282b003..c360bf47a 100644 --- a/include/boost/geometry/policies/relate/intersection_points.hpp +++ b/include/boost/geometry/policies/relate/intersection_points.hpp @@ -16,6 +16,7 @@ #include #include +#include #include #include #include @@ -40,9 +41,6 @@ struct segments_intersection_points S1, S2, CalculationType >::type coordinate_type; - // Get the same type, but at least a double - typedef typename select_most_precise::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 @@ struct segments_intersection_points 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 @@ struct segments_intersection_points // 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(dx1, dy1, dx2, dy2); + promoted_type const da = geometry::determinant(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; diff --git a/include/boost/geometry/strategies/cartesian/area_surveyor.hpp b/include/boost/geometry/strategies/cartesian/area_surveyor.hpp index c9a3e22c4..6f325a7a2 100644 --- a/include/boost/geometry/strategies/cartesian/area_surveyor.hpp +++ b/include/boost/geometry/strategies/cartesian/area_surveyor.hpp @@ -16,9 +16,11 @@ #include -#include -#include +#include +#include +#include +#include namespace boost { namespace geometry @@ -82,7 +84,8 @@ private : inline return_type area() const { return_type result = sum; - result /= 2; + return_type const two = 2; + result /= two; return result; } }; @@ -96,7 +99,7 @@ public : 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(p2, p1); } static inline return_type result(summation const& state) diff --git a/include/boost/geometry/strategies/cartesian/cart_intersect.hpp b/include/boost/geometry/strategies/cartesian/cart_intersect.hpp index d0beeae5c..302055bbd 100644 --- a/include/boost/geometry/strategies/cartesian/cart_intersect.hpp +++ b/include/boost/geometry/strategies/cartesian/cart_intersect.hpp @@ -16,6 +16,7 @@ #include #include +#include #include #include @@ -199,13 +200,19 @@ struct relate_cartesian_segments 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 : " diff --git a/include/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp b/include/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp index e696581cd..5614ec2b1 100644 --- a/include/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp +++ b/include/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp @@ -19,6 +19,7 @@ #include #include +#include #include #include #include @@ -177,7 +178,7 @@ public : calculation_type const y1 = boost::numeric_cast(get<1>(p1)); calculation_type const x2 = boost::numeric_cast(get<0>(p2)); calculation_type const y2 = boost::numeric_cast(get<1>(p2)); - calculation_type const ai = x1 * y2 - x2 * y1; + calculation_type const ai = geometry::cross_product2(p1, p2); state.count++; state.sum_a2 += ai; state.sum_x += ai * (x1 + x2); diff --git a/include/boost/geometry/strategies/cartesian/side_by_triangle.hpp b/include/boost/geometry/strategies/cartesian/side_by_triangle.hpp index 7c4a9d651..a7851d05d 100644 --- a/include/boost/geometry/strategies/cartesian/side_by_triangle.hpp +++ b/include/boost/geometry/strategies/cartesian/side_by_triangle.hpp @@ -17,10 +17,9 @@ #include #include +#include #include - #include - #include @@ -66,7 +65,6 @@ public : 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 @@ public : 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 + ( + dx, dy, + dpx, dpy + ); promoted_type const zero = promoted_type(); return math::equals(s, zero) ? 0