X-Git-Url: https://git.dogcows.com/gitweb?p=chaz%2Fyoink;a=blobdiff_plain;f=src%2FMoof%2Fcml%2Fmathlib%2Fvector_misc.h;fp=src%2FMoof%2Fcml%2Fmathlib%2Fvector_misc.h;h=6553342ba8dbce52d55bfe21a413616a8db0d196;hp=0000000000000000000000000000000000000000;hb=c2321281bf12a7efaedde930422c7ddbc92080d4;hpb=87bc17e55b0c1dc73ecc66df856d3f08fd7a7724 diff --git a/src/Moof/cml/mathlib/vector_misc.h b/src/Moof/cml/mathlib/vector_misc.h new file mode 100644 index 0000000..6553342 --- /dev/null +++ b/src/Moof/cml/mathlib/vector_misc.h @@ -0,0 +1,304 @@ +/* -*- C++ -*- ------------------------------------------------------------ + +Copyright (c) 2007 Jesse Anders and Demian Nave http://cmldev.net/ + +The Configurable Math Library (CML) is distributed under the terms of the +Boost Software License, v1.0 (see cml/LICENSE for details). + + *-----------------------------------------------------------------------*/ +/** @file + * @brief + */ + +#ifndef vector_misc_h +#define vector_misc_h + +#include + +/* Miscellaneous vector functions. */ + +namespace cml { + +/* Function to project a vector v onto a hyperplane specified by a unit-length + * normal n. + * + * @todo: Clean up promotion code. + */ + +template < class VecT_1, class VecT_2 > +typename detail::CrossPromote::promoted_vector +project_to_hplane(const VecT_1& v, const VecT_2& n) +{ + typedef typename detail::CrossPromote::promoted_vector + result_type; + + result_type result; + et::detail::Resize(result, v.size()); + + result = v - dot(v,n) * n; + return result; +} + +/* Return a vector perpendicular (CCW) to a 2D vector. */ +template < class VecT > vector< typename VecT::value_type, fixed<2> > +perp(const VecT& v) +{ + typedef vector< typename VecT::value_type, fixed<2> > temporary_type; + + /* Checking */ + detail::CheckVec2(v); + + return temporary_type(-v[1],v[0]); +} + +/* @todo: unit_cross() and cross_cardinal() should probably go in + * vector_products.h, but I'm trying to avoid modifying the existing codebase + * for now. + */ + +/** Return normalized cross product of two vectors */ +template< class LeftT, class RightT > +typename detail::CrossPromote::promoted_vector +unit_cross(const LeftT& left, const RightT& right) { + /* @todo: This will probably break with dynamic<> vectors */ + return normalize(cross(left,right)); +} + +/** Return the cross product of v and the i'th cardinal basis vector */ +template < class VecT > vector< typename VecT::value_type, fixed<3> > +cross_cardinal(const VecT& v, size_t i) +{ + typedef vector< typename VecT::value_type, fixed<3> > vector_type; + typedef typename vector_type::value_type value_type; + + /* Checking */ + detail::CheckVec3(v); + detail::CheckIndex3(i); + + size_t j, k; + cyclic_permutation(i, i, j, k); + vector_type result; + result[i] = value_type(0); + result[j] = v[k]; + result[k] = -v[j]; + return result; +} + +/** Return the cross product of the i'th cardinal basis vector and v */ +template < class VecT > vector< typename VecT::value_type, fixed<3> > +cross_cardinal(size_t i, const VecT& v) +{ + typedef vector< typename VecT::value_type, fixed<3> > vector_type; + typedef typename vector_type::value_type value_type; + + /* Checking */ + detail::CheckVec3(v); + detail::CheckIndex3(i); + + size_t j, k; + cyclic_permutation(i, i, j, k); + vector_type result; + result[i] = value_type(0); + result[j] = -v[k]; + result[k] = v[j]; + return result; +} + +/** Rotate a 3D vector v about a unit-length vector n. */ +template< class VecT_1, class VecT_2, typename Real > +vector< + typename et::ScalarPromote< + typename VecT_1::value_type, + typename VecT_2::value_type + >::type, + fixed<3> +> +rotate_vector(const VecT_1& v, const VecT_2& n, Real angle) +{ + typedef vector< + typename et::ScalarPromote< + typename VecT_1::value_type, + typename VecT_2::value_type + >::type, + fixed<3> + > result_type; + + /* Checking */ + detail::CheckVec3(v); + detail::CheckVec3(n); + + result_type parallel = dot(v,n)*n; + return ( + std::cos(angle)*(v-parallel) + std::sin(angle)*cross(n,v) + parallel + ); +} + +/** Rotate a 2D vector v about a unit-length vector n. */ +template< class VecT, typename Real > +vector< typename VecT::value_type, fixed<2> > +rotate_vector_2D(const VecT& v, Real angle) +{ + typedef vector< typename VecT::value_type, fixed<2> > result_type; + typedef typename result_type::value_type value_type; + + /* Checking */ + detail::CheckVec2(v); + + value_type s = std::sin(angle); + value_type c = std::cos(angle); + + return result_type(c * v[0] - s * v[1], s * v[0] + c * v[1]); +} + +/** Random unit 3D or 2D vector + * + * @todo: This is just placeholder code for what will be a more thorough + * 'random unit' implementation: + * + * - All dimensions will be handled uniformly if practical, perhaps through + * a normal distrubution PRNG. + * + * - Failing that (or perhaps even in this case), dimensions 2 and 3 will be + * dispatched to special-case code, most likely implementing the algorithms + * below. + * + * - Like the utility random functions, the option of using one's own PRGN + * will be made available. + * + * @todo: Once N-d random vectors are supported, add a 'random unit + * quaternion' function that wraps a call to random_unit() with a 4D vector as + * the argument. + */ +template < typename E, class A > void +random_unit(vector& v) +{ + typedef vector vector_type; + typedef typename vector_type::value_type value_type; + + switch (v.size()) { + case 3: + { + vector< E, fixed<3> > temp; + spherical_to_cartesian( + value_type(1), + value_type(random_unit() * constants::two_pi()), + acos_safe(random_real(value_type(-1),value_type(1))), + 2, + colatitude, + temp + ); + v[0] = temp[0]; + v[1] = temp[1]; + v[2] = temp[2]; + break; + } + case 2: + { + vector< E, fixed<2> > temp; + polar_to_cartesian( + value_type(1), + value_type(random_unit() * constants::two_pi()), + temp + ); + v[0] = temp[0]; + v[1] = temp[1]; + break; + } + default: + throw std::invalid_argument( + "random_unit() for N-d vectors not implemented yet"); + break; + } +} + +/* Random vector within a given angle of a unit-length axis, i.e. in a cone + * (3D) or wedge (2D). + * + * The same notes the appear above apply here too, more or less. One + * difference is that this is really only useful in 2D and 3D (presumably), so + * we'll probably just do a compile- or run-time dispatch as appropriate. + * + * Also, there may be a better algorithm for generating a random unit vector + * in a cone; need to look into that. + * + * All of this 'temp' stuff is because there's no compile-time dispatch for + * 3D and 2D vectors, but that'll be fixed soon. + */ + +template < typename E, class A, class VecT > void +random_unit(vector& v, const VecT& axis, E theta) +{ + typedef vector vector_type; + typedef typename vector_type::value_type value_type; + + switch (v.size()) { + case 3: + { + vector< E, fixed<3> > temp, n, temp_axis; + temp_axis[0] = axis[0]; + temp_axis[1] = axis[1]; + temp_axis[2] = axis[2]; + + /* @todo: Function for finding 'any perpendicular vector'? */ + n = axis_3D(cml::index_of_min_abs(axis[0],axis[1],axis[2])); + n = cross(n,temp_axis); + + /* Rotate v 'away from' the axis by a random angle in the range + * [-theta,theta] + */ + temp = rotate_vector(temp_axis,n,random_real(-theta,theta)); + + /* Rotate v about the axis by a random angle in the range [-pi,pi] + */ + temp = rotate_vector( + temp, + temp_axis, + random_real( + -constants::pi(), + constants::pi() + ) + ); + + v[0] = temp[0]; + v[1] = temp[1]; + v[2] = temp[2]; + break; + } + case 2: + { + vector< E, fixed<2> > temp, temp_axis; + temp_axis[0] = axis[0]; + temp_axis[1] = axis[1]; + temp = rotate_vector_2D(temp_axis, random_real(-theta,theta)); + v[0] = temp[0]; + v[1] = temp[1]; + break; + } + default: + throw std::invalid_argument( + "random_unit(v,axis,theta) only implemented for 2D and 3D"); + break; + } +} + +/* NEW: Manhattan distance */ + +template< class VecT_1, class VecT_2 > +typename detail::DotPromote< VecT_1, VecT_2 >::promoted_scalar +manhattan_distance(const VecT_1& v1, const VecT_2& v2) { + /* Check that a promotion exists */ + typedef typename et::VectorPromote< + VecT_1,VecT_2>::temporary_type promoted_vector; + + typedef typename detail::DotPromote< VecT_1, VecT_2 >::promoted_scalar scalar_type; + + scalar_type sum = scalar_type(0); + for (size_t i = 0; i < v1.size(); ++i) { + sum += std::fabs(v2[i]-v1[i]); + } + return sum; +} + +} // namespace cml + +#endif