boost/geometry
2
0
Fork 0
mirror of https://github.com/boostorg/geometry.git synced 2026-02-12 00:02:09 +00:00
Code Issues Projects Releases Wiki Activity

Introduced cross-product for area,centroid,side,intersection(determinant,direction,relation)

Browse Source 
[SVN r76755]
This commit is contained in:
Barend Gehrels
2012-01-28 18:29:47 +00:00
parent 10b649c234
commit 4594a65da4
7 changed files with 137 additions and 40 deletions
Download Patch File Download Diff File Expand all files Collapse all files

89
include/boost/geometry/arithmetic/cross_product.hpp
View File

@@ -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 @@ struct cross_product<P1, P2, 3>
}
};
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 @@ struct cross_product<P1, P2, 3>
// -- 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 @@ 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 <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

30
include/boost/geometry/policies/relate/direction.hpp
View File

@@ -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 @@ 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<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 @@ private :
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 @@ private :
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 @@ 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);
}

20
include/boost/geometry/policies/relate/intersection_points.hpp
View File

@@ -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 @@ struct segments_intersection_points
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 @@ 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<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;

11
include/boost/geometry/strategies/cartesian/area_surveyor.hpp
View File

@@ -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 @@ 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<return_type>(p2, p1);
}
static inline return_type result(summation const& state)

13
include/boost/geometry/strategies/cartesian/cart_intersect.hpp
View File

@@ -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 @@ 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 : "

3
include/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp
View File

@@ -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 @@ public :
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);

11
include/boost/geometry/strategies/cartesian/side_by_triangle.hpp
View File

@@ -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 @@ 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<promoted_type>
(
dx, dy,
dpx, dpy
);
promoted_type const zero = promoted_type();
return math::equals(s, zero) ? 0