*/
template<typename Float>
struct constants {
- static Float pi() { return M_PI; }
- static Float two_pi() { return 2.*M_PI; }
- static Float inv_pi() { return 1./M_PI; }
- static Float inv_two_pi() { return 1./(2.*M_PI); }
- static Float pi_over_2() { return M_PI/2.; }
- static Float pi_over_4() { return M_PI/4.; }
- static Float deg_per_rad() { return 180./M_PI; }
- static Float rad_per_deg() { return M_PI/180.; }
- static Float sqrt_2() { return M_SQRT2; }
- static Float sqrt_3() { return 1.732050807568877293527446341505; }
- static Float sqrt_5() { return 2.236067977499789696409173668731; }
- static Float sqrt_6() { return 2.449489742783178098197284074705; }
- static Float e() { return M_E; }
+ static Float pi() { return Float(M_PI); }
+ static Float two_pi() { return Float(2.*M_PI); }
+ static Float inv_pi() { return Float(1./M_PI); }
+ static Float inv_two_pi() { return Float(1./(2.*M_PI)); }
+ static Float pi_over_2() { return Float(M_PI/2.); }
+ static Float pi_over_4() { return Float(M_PI/4.); }
+ static Float deg_per_rad() { return Float(180./M_PI); }
+ static Float rad_per_deg() { return Float(M_PI/180.); }
+
+ static Float sqrt_2() { return Float(M_SQRT2); }
+ static Float sqrt_3() { return Float(1.732050807568877293527446341505); }
+ static Float sqrt_5() { return Float(2.236067977499789696409173668731); }
+ static Float sqrt_6() { return Float(2.449489742783178098197284074705); }
+
+ static Float e() { return Float(M_E); }
};
#else
enum { is_true = true, is_false = false };
};
+/** Remove a reference qualifier from a type. */
+template<typename T> struct remove_reference {
+ template<typename Q, typename Dummy> struct helper {
+ typedef Q type;
+ };
+
+ template<typename Q> struct helper<Q&, void> {
+ typedef Q type;
+ };
+
+ template<typename Q> struct helper<const Q&, void> {
+ typedef const Q type;
+ };
+
+ typedef typename helper<T,void>::type type;
+};
+
+/** Remove a const qualifier from a type. */
+template<typename T> struct remove_const {
+ template<typename Q, typename Dummy> struct helper {
+ typedef Q type;
+ };
+
+ template<typename Q> struct helper<const Q, void> {
+ typedef Q type;
+ };
+
+ typedef typename helper<T,void>::type type;
+};
+
} // namespace cml
#endif
#ifndef scalar_promotions_h
#define scalar_promotions_h
+#include <complex>
#include <cml/core/cml_meta.h>
namespace cml {
namespace et {
+// #define CML_USE_OLD_SCALAR_PROMOTIONS
+#if !defined(CML_USE_OLD_SCALAR_PROMOTIONS)
+
+/* The type promotion code below is a slightly modified version of:
+ * http://ubiety.uwaterloo.ca/~tveldhui/papers/techniques/techniques01.html
+ */
+namespace detail {
+
+template<class T>
+struct precision_trait {
+ enum { precisionRank = 0,
+ knowPrecisionRank = 0 };
+};
+
+#define DECLARE_PRECISION(T,rank) \
+ template<> \
+ struct precision_trait< T > { \
+ enum { precisionRank = rank, \
+ knowPrecisionRank = 1 }; \
+ };
+
+DECLARE_PRECISION(int,100)
+DECLARE_PRECISION(unsigned int,200)
+DECLARE_PRECISION(long,300)
+DECLARE_PRECISION(unsigned long,400)
+
+DECLARE_PRECISION(long long,425)
+DECLARE_PRECISION(unsigned long long,475)
+
+DECLARE_PRECISION(float,500)
+DECLARE_PRECISION(double,600)
+DECLARE_PRECISION(long double,700)
+DECLARE_PRECISION(std::complex<float>,800)
+DECLARE_PRECISION(std::complex<double>,900)
+DECLARE_PRECISION(std::complex<long double>,1000)
+
+template<class T>
+struct autopromote_trait {
+ typedef T T_numtype;
+};
+
+#define DECLARE_AUTOPROMOTE(T1,T2) \
+ template<> \
+ struct autopromote_trait<T1> { \
+ typedef T2 T_numtype; \
+ };
+
+// These are the odd cases where small integer types
+// are automatically promoted to int or unsigned int for
+// arithmetic.
+DECLARE_AUTOPROMOTE(bool, int)
+DECLARE_AUTOPROMOTE(char, int)
+DECLARE_AUTOPROMOTE(unsigned char, int)
+DECLARE_AUTOPROMOTE(short int, int)
+DECLARE_AUTOPROMOTE(short unsigned int, unsigned int)
+
+template<class T1, class T2, int promoteToT1>
+struct promote2 {
+ typedef T1 T_promote;
+};
+
+template<class T1, class T2>
+struct promote2<T1,T2,0> {
+ typedef T2 T_promote;
+};
+
+template<class T1_orig, class T2_orig>
+struct promote_trait {
+ // Handle promotion of small integers to int/unsigned int
+ typedef typename autopromote_trait<T1_orig>::T_numtype T1;
+ typedef typename autopromote_trait<T2_orig>::T_numtype T2;
+
+ // True if T1 is higher ranked
+ enum {
+ T1IsBetter =
+ (int) precision_trait<T1>::precisionRank >
+ (int) precision_trait<T2>::precisionRank,
+
+ // True if we know ranks for both T1 and T2
+ knowBothRanks =
+ precision_trait<T1>::knowPrecisionRank
+ && precision_trait<T2>::knowPrecisionRank,
+
+ // True if we know T1 but not T2
+ knowT1butNotT2 = precision_trait<T1>::knowPrecisionRank
+ && !(precision_trait<T2>::knowPrecisionRank),
+
+ // True if we know T2 but not T1
+ knowT2butNotT1 = precision_trait<T2>::knowPrecisionRank
+ && !(precision_trait<T1>::knowPrecisionRank),
+
+ // True if T1 is bigger than T2
+ T1IsLarger = sizeof(T1) >= sizeof(T2),
+
+ // We know T1 but not T2: true
+ // We know T2 but not T1: false
+ // Otherwise, if T1 is bigger than T2: true
+ defaultPromotion = knowT1butNotT2 ? false :
+ (knowT2butNotT1 ? true : T1IsLarger)
+ };
+
+ // If we have both ranks, then use them.
+ // If we have only one rank, then use the unknown type.
+ // If we have neither rank, then promote to the larger type.
+
+ enum {
+ promoteToT1 = (knowBothRanks ? T1IsBetter : defaultPromotion)
+ ? 1 : 0
+ };
+
+ typedef typename promote2<T1,T2,promoteToT1>::T_promote T_promote;
+};
+
+} // namespace detail
+
+/** Defers to detail::promote_trait<>. */
+template<class E1, class E2> struct ScalarPromote
+{
+ typedef typename detail::promote_trait<E1,E2>::T_promote type;
+};
+
+#else
+
namespace detail {
/** @class IntPromote
*/
template<class E1_in, class E2_in> struct ScalarPromote
{
-
/* Integral-promote the types (if possible). */
typedef typename detail::IntPromote<E1_in>::result E1;
typedef typename detail::IntPromote<E2_in>::result E2;
same_type<float_filter,void>::is_true,
int_filter, float_filter>::result type;
};
+#endif
} // namespace et
} // namespace cml
planes[5][3] = m.basis_element(3,3) - m.basis_element(3,2);
/* Near: [03+02, 13+12, 23+22, 33+32] : [02, 12, 22, 32] */
-
- if (z_clip == z_clip_neg_one) {
- planes[4][0] = m.basis_element(0,3) + m.basis_element(0,2);
- planes[4][1] = m.basis_element(1,3) + m.basis_element(1,2);
- planes[4][2] = m.basis_element(2,3) + m.basis_element(2,2);
- planes[4][3] = m.basis_element(3,3) + m.basis_element(3,2);
- } else { // z_clip == z_clip_zero
- planes[4][0] = m.basis_element(0,2);
- planes[4][1] = m.basis_element(1,2);
- planes[4][2] = m.basis_element(2,2);
- planes[4][3] = m.basis_element(3,2);
- }
+ extract_near_frustum_plane(m, planes[4], z_clip);
/* @todo: This will be handled by the plane class */
if (normalize) {
}
}
+/** Extract the near plane of a frustum given a concatenated modelview and
+ * projection matrix with the given near z-clipping range. The plane is
+ * not normalized.
+ *
+ * @note The plane is in ax+by+cz+d = 0 form.
+ *
+ * @warning The matrix is assumed to be a homogeneous transformation
+ * matrix.
+ */
+template < class MatT, class PlaneT > void
+extract_near_frustum_plane(
+ const MatT& m,
+ PlaneT& plane,
+ ZClip z_clip
+ )
+{
+ /* Near: [03+02, 13+12, 23+22, 33+32] : [02, 12, 22, 32] */
+ if (z_clip == z_clip_neg_one) {
+ plane[0] = m.basis_element(0,3) + m.basis_element(0,2);
+ plane[1] = m.basis_element(1,3) + m.basis_element(1,2);
+ plane[2] = m.basis_element(2,3) + m.basis_element(2,2);
+ plane[3] = m.basis_element(3,3) + m.basis_element(3,2);
+ } else { // z_clip == z_clip_zero
+ plane[0] = m.basis_element(0,2);
+ plane[1] = m.basis_element(1,2);
+ plane[2] = m.basis_element(2,2);
+ plane[3] = m.basis_element(3,2);
+ }
+}
+
namespace detail {
/* This is currently only in support of finding the corners of a frustum.
* euler_order_zyx
* euler_order_zxz
* euler_order_zyz
+ *
+ * e.g. euler_order_xyz means compute the column-basis rotation matrix
+ * equivalent to R_x * R_y * R_z, where R_i is the rotation matrix above
+ * axis i (the row-basis matrix would be R_z * R_y * R_x).
*/
-
template < typename E, class A, class B, class L > void
matrix_rotation_euler(matrix<E,A,B,L>& m, E angle_0, E angle_1, E angle_2,
EulerOrder order)
}
}
+/** Build a matrix of derivatives of Euler angles about the specified axis.
+ *
+ * The rotation derivatives are applied about the cardinal axes in the
+ * order specified by the 'order' argument, where 'order' is one of the
+ * following enumerants:
+ *
+ * euler_order_xyz
+ * euler_order_xzy
+ * euler_order_yzx
+ * euler_order_yxz
+ * euler_order_zxy
+ * euler_order_zyx
+ *
+ * e.g. euler_order_xyz means compute the column-basis rotation matrix
+ * equivalent to R_x * R_y * R_z, where R_i is the rotation matrix above
+ * axis i (the row-basis matrix would be R_z * R_y * R_x).
+ *
+ * The derivative is taken with respect to the specified 'axis', which is
+ * the position of the axis in the triple; e.g. if order = euler_order_xyz,
+ * then axis = 0 would mean take the derivative with respect to x. Note
+ * that repeated axes are not currently supported.
+ */
+template < typename E, class A, class B, class L > void
+matrix_rotation_euler_derivatives(
+ matrix<E,A,B,L>& m, int axis, E angle_0, E angle_1, E angle_2,
+ EulerOrder order)
+{
+ typedef matrix<E,A,B,L> matrix_type;
+ typedef typename matrix_type::value_type value_type;
+
+ /* Checking */
+ detail::CheckMatLinear3D(m);
+
+ identity_transform(m);
+
+ size_t i, j, k;
+ bool odd, repeat;
+ detail::unpack_euler_order(order, i, j, k, odd, repeat);
+ if(repeat) throw std::invalid_argument(
+ "matrix_rotation_euler_derivatives does not support repeated axes");
+
+ if (odd) {
+ angle_0 = -angle_0;
+ angle_1 = -angle_1;
+ angle_2 = -angle_2;
+ }
+
+ value_type s0 = std::sin(angle_0);
+ value_type c0 = std::cos(angle_0);
+ value_type s1 = std::sin(angle_1);
+ value_type c1 = std::cos(angle_1);
+ value_type s2 = std::sin(angle_2);
+ value_type c2 = std::cos(angle_2);
+
+ value_type s0s2 = s0 * s2;
+ value_type s0c2 = s0 * c2;
+ value_type c0s2 = c0 * s2;
+ value_type c0c2 = c0 * c2;
+
+ if(axis == 0) {
+ m.set_basis_element(i,i, 0. );
+ m.set_basis_element(i,j, 0. );
+ m.set_basis_element(i,k, 0. );
+ m.set_basis_element(j,i, s1 * c0*c2 + s0*s2);
+ m.set_basis_element(j,j, s1 * c0*s2 - s0*c2);
+ m.set_basis_element(j,k, c0 * c1 );
+ m.set_basis_element(k,i,-s1 * s0*c2 + c0*s2);
+ m.set_basis_element(k,j,-s1 * s0*s2 - c0*c2);
+ m.set_basis_element(k,k,-s0 * c1 );
+ } else if(axis == 1) {
+ m.set_basis_element(i,i,-s1 * c2 );
+ m.set_basis_element(i,j,-s1 * s2 );
+ m.set_basis_element(i,k,-c1 );
+ m.set_basis_element(j,i, c1 * s0*c2 );
+ m.set_basis_element(j,j, c1 * s0*s2 );
+ m.set_basis_element(j,k,-s0 * s1 );
+ m.set_basis_element(k,i, c1 * c0*c2 );
+ m.set_basis_element(k,j, c1 * c0*s2 );
+ m.set_basis_element(k,k,-c0 * s1 );
+ } else if(axis == 2) {
+ m.set_basis_element(i,i,-c1 * s2 );
+ m.set_basis_element(i,j, c1 * c2 );
+ m.set_basis_element(i,k, 0. );
+ m.set_basis_element(j,i,-s1 * s0*s2 - c0*c2);
+ m.set_basis_element(j,j, s1 * s0*c2 - c0*s2);
+ m.set_basis_element(j,k, 0. );
+ m.set_basis_element(k,i,-s1 * c0*s2 + s0*c2);
+ m.set_basis_element(k,j, s1 * c0*c2 + s0*s2);
+ m.set_basis_element(k,k, 0. );
+ }
+}
+
//////////////////////////////////////////////////////////////////////////////
// 3D rotation to align with a vector, multiple vectors, or the view plane
//////////////////////////////////////////////////////////////////////////////
);
}
+/** Get the translation of a 3D affine transform */
+template < class MatT > void
+matrix_get_translation(
+ const MatT& m,
+ typename MatT::value_type& t1,
+ typename MatT::value_type& t2,
+ typename MatT::value_type& t3
+ )
+{
+ typedef typename MatT::value_type value_type;
+ typedef vector< value_type, fixed<3> > vector_type;
+
+ /* Checking */
+ detail::CheckMatAffine3D(m);
+
+ t1 = m.basis_element(3,0);
+ t2 = m.basis_element(3,1);
+ t3 = m.basis_element(3,2);
+}
+
/** Get the translation of a 2D affine transform */
template < class MatT > vector< typename MatT::value_type, fixed<2> >
matrix_get_translation_2D(const MatT& m)
return vector_type(m.basis_element(2,0), m.basis_element(2,1));
}
+/** Get the translation of a 2D affine transform */
+template < class MatT > void
+matrix_get_translation_2D(
+ const MatT& m,
+ typename MatT::value_type& t1,
+ typename MatT::value_type& t2
+ )
+{
+ typedef typename MatT::value_type value_type;
+ typedef vector< value_type, fixed<2> > vector_type;
+
+ /* Checking */
+ detail::CheckMatAffine2D(m);
+
+ t1 = m.basis_element(2,0);
+ t2 = m.basis_element(2,1);
+}
+
//////////////////////////////////////////////////////////////////////////////
// Function for getting the translation of a 3D view matrix
//////////////////////////////////////////////////////////////////////////////
typedef vector< double, fixed<4> > vector4d;
/* fixed-size matrices */
+
+typedef matrix< int, fixed<2,2> > matrix22i;
+typedef matrix< float, fixed<2,2> > matrix22f;
+typedef matrix< double, fixed<2,2> > matrix22d;
+
typedef matrix< int, fixed<2,2>, row_basis, row_major > matrix22i_r;
typedef matrix< int, fixed<2,2>, col_basis, col_major > matrix22i_c;
typedef matrix< float, fixed<2,2>, row_basis, row_major > matrix22f_r;
typedef matrix< double, fixed<2,2>, row_basis, row_major > matrix22d_r;
typedef matrix< double, fixed<2,2>, col_basis, col_major > matrix22d_c;
+
+typedef matrix< int, fixed<3,3> > matrix33i;
+typedef matrix< float, fixed<3,3> > matrix33f;
+typedef matrix< double, fixed<3,3> > matrix33d;
+
typedef matrix< int, fixed<3,3>, row_basis, row_major > matrix33i_r;
typedef matrix< int, fixed<3,3>, col_basis, col_major > matrix33i_c;
typedef matrix< float, fixed<3,3>, row_basis, row_major > matrix33f_r;
typedef matrix< double, fixed<3,3>, row_basis, row_major > matrix33d_r;
typedef matrix< double, fixed<3,3>, col_basis, col_major > matrix33d_c;
+
+typedef matrix< int, fixed<4,4> > matrix44i;
+typedef matrix< float, fixed<4,4> > matrix44f;
+typedef matrix< double, fixed<4,4> > matrix44d;
+
typedef matrix< int, fixed<4,4>, row_basis, row_major > matrix44i_r;
typedef matrix< int, fixed<4,4>, col_basis, col_major > matrix44i_c;
typedef matrix< float, fixed<4,4>, row_basis, row_major > matrix44f_r;
typedef matrix< double, fixed<4,4>, row_basis, row_major > matrix44d_r;
typedef matrix< double, fixed<4,4>, col_basis, col_major > matrix44d_c;
+
typedef matrix< int, fixed<3,2>, row_basis, row_major > matrix32i_r;
typedef matrix< float, fixed<3,2>, row_basis, row_major > matrix32f_r;
typedef matrix< double, fixed<3,2>, row_basis, row_major > matrix32d_r;
typedef matrix< float, fixed<3,4>, col_basis, col_major > matrix34f_c;
typedef matrix< double, fixed<3,4>, col_basis, col_major > matrix34d_c;
+
/* quaternions */
typedef quaternion<float, fixed<>,vector_first,negative_cross>
quaternionf_n;
quaterniond_n;
typedef quaternion<double,fixed<>,vector_first,positive_cross>
quaterniond_p;
+typedef quaternion<float> quaternionf;
+typedef quaternion<double> quaterniond;
+
/* dynamically resizable vectors */
typedef vector< int, dynamic<> > vectori;
typedef vector< float, dynamic<> > vectorf;
typedef vector< double, dynamic<> > vectord;
+
/* dynamically resizable matrices */
+typedef matrix< int, dynamic<> > matrixi;
+typedef matrix< float, dynamic<> > matrixf;
+typedef matrix< double, dynamic<> > matrixd;
+
typedef matrix< int, dynamic<>, row_basis, row_major > matrixi_r;
typedef matrix< int, dynamic<>, col_basis, col_major > matrixi_c;
typedef matrix< float, dynamic<>, row_basis, row_major > matrixf_r;
typedef matrix< double, dynamic<>, row_basis, row_major > matrixd_r;
typedef matrix< double, dynamic<>, col_basis, col_major > matrixd_c;
+
/* constants */
typedef constants<float> constantsf;
typedef constants<double> constantsd;
*-----------------------------------------------------------------------*/
/** @file
- * @brief
+ *
+ * Functions for orthonormalizing a set of basis vectors in 3D or 2D, and for
+ * constructing an orthonormal basis given various input parameters.
*/
#ifndef vector_ortho_h
#include <cml/mathlib/vector_misc.h>
#include <cml/mathlib/misc.h>
-/* Functions for orthonormalizing a set of basis vector in 3D or 2D, and for
- * constructing an orthonormal basis given various input parameters.
- */
-
namespace cml {
-/* Orthonormalize 3 basis vectors in R3.
+//////////////////////////////////////////////////////////////////////////////
+// Orthonormalization in 3D and 2D
+//////////////////////////////////////////////////////////////////////////////
+
+
+/** Orthonormalize 3 basis vectors in R3.
*
* Called with the default values, this function performs a single Gram-
* Schmidt step to orthonormalize the input vectors. By default, the direction
* In most cases, the default arguments can be ignored, leaving only the three
* input vectors.
*/
-
-//////////////////////////////////////////////////////////////////////////////
-// Orthonormalization in 3D and 2D
-//////////////////////////////////////////////////////////////////////////////
-
template < typename E, class A > void
orthonormalize(vector<E,A>& v0, vector<E,A>& v1, vector<E,A>& v2,
size_t stable_axis = 2, size_t num_iter = 0, E s = E(1))
v2 = v[2];
}
-/* Orthonormalize 2 basis vectors in R2 */
+/** Orthonormalize 2 basis vectors in R2 */
template < typename E, class A > void
orthonormalize(vector<E,A>& v0, vector<E,A>& v1,
size_t stable_axis = 0, size_t num_iter = 0, E s = E(1))
// Orthonormal basis construction in 3D and 2D
//////////////////////////////////////////////////////////////////////////////
-/* This version of orthonormal_basis() ultimately does the work for all
+/** This version of orthonormal_basis() ultimately does the work for all
* orthonormal_basis_*() functions. Given input vectors 'align' and
- * 'reference', and an order 'axis_order_<i><j><k>', it constructs an
+ * 'reference', and an order 'axis_order_\<i\>\<j\>\<k\>', it constructs an
* orthonormal basis such that the i'th basis vector is aligned with (parallel
* to and pointing in the same direction as) 'align', and the j'th basis
* vector is maximally aligned with 'reference'. The k'th basis vector is
* chosen such that the basis has a determinant of +1.
*
- * Note that the algorithm fails when 'align' is nearly parallel to
+ * @note The algorithm fails when 'align' is nearly parallel to
* 'reference'; this should be checked for and handled externally if it's a
* case that may occur.
- */
-
-/* Note: This is an example of the 'non-const argument modification
+ *
+ * @internal This is an example of the 'non-const argument modification
* invalidates expression' gotcha. If x, y or z were to be assigned to before
* we were 'done' with align and reference, and if one of them were the same
* object as align or reference, then the algorithm could fail. As is the
* basis vectors are assigned at the end of the function from a temporary
* array, so all is well.
*/
-
template < class VecT_1, class VecT_2, typename E, class A > void
orthonormal_basis(
const VecT_1& align,
z = axis[2];
}
-/* This version of orthonormal_basis() constructs in arbitrary basis given a
+/** This version of orthonormal_basis() constructs in arbitrary basis given a
* vector with which to align the i'th basis vector. To avoid the failure
* case, the reference vector is always chosen so as to not be parallel to
* 'align'. This means the algorithm will always generate a valid basis, which
* basis will likely 'pop' as the alignment vector changes, and so may not be
* suitable for billboarding or other similar applications.
*/
-
template < class VecT, typename E, class A >
void orthonormal_basis(
const VecT& align,
);
}
-/* orthonormal_basis_axial() generates a basis in which the j'th basis vector
+/** orthonormal_basis_axial() generates a basis in which the j'th basis vector
* is aligned with 'axis' and the i'th basis vector is maximally aligned (as
* 'aligned as possible') with 'align'. This can be used for e.g. axial
* billboarding for, say, trees or beam effects.
* of orthonormal_basis(), with the parameters adjusted so that the alignment
* is axial.
*
- * With this algorithm the failure case is when 'align' and 'axis' are nearly
- * parallel; if this is likely, it should be checked for and handled
- * externally.
+ * @note With this algorithm the failure case is when 'align' and 'axis'
+ * are nearly parallel; if this is likely, it should be checked for and
+ * handled externally.
*/
-
template < class VecT_1, class VecT_2, typename E, class A >
void orthonormal_basis_axial(
const VecT_1& align,
detail::swap_axis_order(order));
}
-/* orthonormal_basis_viewplane() builds a basis aligned with a viewplane, as
+/** orthonormal_basis_viewplane() builds a basis aligned with a viewplane, as
* extracted from the input view matrix. The function takes into account the
* handedness of the input view matrix and orients the basis accordingly.
*
- * The generated basis will always be valid.
+ * @note The generated basis will always be valid.
*/
template < class MatT, typename E, class A >
void orthonormal_basis_viewplane(
view_matrix,x,y,z,right_handed,order);
}
-/* Build a 2D orthonormal basis. */
+/** Build a 2D orthonormal basis. */
template < class VecT, typename E, class A >
void orthonormal_basis_2D(
const VecT& align,
set_basis_element(i,j,s,basis_orient());
}
+ /** Set the matrix row from the given vector. */
+ void set_row(size_t i, const row_vector_type& row) {
+ for(size_t j = 0; j < this->cols(); ++ j) (*this)(i,j) = row[j];
+ }
+
+ /** Set the matrix column from the given vector. */
+ void set_col(size_t j, const col_vector_type& col) {
+ for(size_t i = 0; i < this->rows(); ++ i) (*this)(i,j) = col[i];
+ }
+
public:
set_basis_element(i,j,s,basis_orient());
}
+ /** Set the matrix row from the given vector. */
+ void set_row(size_t i, const row_vector_type& row) {
+ for(size_t j = 0; j < this->cols(); ++ j) (*this)(i,j) = row[j];
+ }
+
+ /** Set the matrix column from the given vector. */
+ void set_col(size_t j, const col_vector_type& col) {
+ for(size_t i = 0; i < this->rows(); ++ i) (*this)(i,j) = col[i];
+ }
+
public:
* responsible for doing any necessary memory management.
*
* @param ptr specify the external pointer.
+ * @param rows the number of rows in the C array.
+ * @param cols the number of columns in the C array.
*
* @throws same as the ArrayType constructor.
*/
set_basis_element(i,j,s,basis_orient());
}
+ /** Set the matrix row from the given vector. */
+ void set_row(size_t i, const row_vector_type& row) {
+ for(size_t j = 0; j < this->cols(); ++ j) (*this)(i,j) = row[j];
+ }
+
+ /** Set the matrix column from the given vector. */
+ void set_col(size_t j, const col_vector_type& col) {
+ for(size_t i = 0; i < this->rows(); ++ i) (*this)(i,j) = col[i];
+ }
+
public:
namespace cml {
namespace et {
-/* Default matrix type promotion template. */
-template<typename LeftT, typename RightT> struct MatrixPromote;
-
-/** Type promotion for two matrix types.
- *
- * @note This always uses the basis orientation of the left-hand matrix.
- * @bug This always uses the basis orientation of the left-hand matrix,
- * which is not always correct.
- */
-template<typename E1, class AT1, typename L1, typename BO1,
- typename E2, class AT2, typename L2, typename BO2>
-struct MatrixPromote< cml::matrix<E1,AT1,BO1,L1>, cml::matrix<E2,AT2,BO2,L2> >
+/** Promote two types to a matrixt type. */
+template<typename LeftT, typename RightT> struct MatrixPromote
{
- /* Promote the arrays: */
- typedef typename ArrayPromote<
- typename cml::matrix<E1,AT1,BO1,L1>::array_type,
+ /* Default matrix type promotion template. */
+ template<typename M1, typename M2> struct MatrixPromoteHelper;
+
+ /** Type promotion for two matrix types.
+ *
+ * @note This always uses the basis orientation of the left-hand matrix.
+ * @bug This always uses the basis orientation of the left-hand matrix,
+ * which is not always correct.
+ */
+ template<typename E1, class AT1, typename L1, typename BO1,
+ typename E2, class AT2, typename L2, typename BO2>
+ struct MatrixPromoteHelper<
+ cml::matrix<E1,AT1,BO1,L1>, cml::matrix<E2,AT2,BO2,L2>
+ >
+ {
+ /* Promote the arrays: */
+ typedef typename ArrayPromote<
+ typename cml::matrix<E1,AT1,BO1,L1>::array_type,
typename cml::matrix<E2,AT2,BO2,L2>::array_type
- >::type promoted_array;
+ >::type promoted_array;
+
+ /* The deduced matrix result type: */
+ typedef cml::matrix<
+ typename promoted_array::value_type,
+ typename promoted_array::generator_type,
+ BO1,
+ typename promoted_array::layout
+ > type;
+
+ /* The deduced temporary type: */
+ typedef typename type::temporary_type temporary_type;
+ };
+
+ /** Type promotion for a matrix and a scalar. */
+ template<typename E, class AT, typename BO, typename L, typename S>
+ struct MatrixPromoteHelper<cml::matrix<E,AT,BO,L>, S>
+ {
+ /* The deduced matrix result type (the array type is the same): */
+ typedef cml::matrix<typename ScalarPromote<E,S>::type, AT, BO, L> type;
+
+ /* The deduced temporary type: */
+ typedef typename type::temporary_type temporary_type;
+ };
+
+ /** Type promotion for a scalar and a matrix. */
+ template<typename S, typename E, class AT, typename BO, typename L>
+ struct MatrixPromoteHelper<S, cml::matrix<E,AT,BO,L> >
+ {
+ /* The deduced matrix result type (the array type is the same): */
+ typedef cml::matrix<typename ScalarPromote<S,E>::type, AT, BO, L> type;
- /* The deduced matrix result type: */
- typedef cml::matrix<
- typename promoted_array::value_type,
- typename promoted_array::generator_type,
- BO1,
- typename promoted_array::layout
- > type;
+ /* The deduced temporary type: */
+ typedef typename type::temporary_type temporary_type;
+ };
+
+ /** Type promotion for outer product. */
+ template<typename E1, class AT1, typename E2, class AT2>
+ struct MatrixPromoteHelper< cml::vector<E1,AT1>, cml::vector<E2,AT2> >
+ {
+ typedef cml::vector<E1,AT1> left_type;
+ typedef cml::vector<E2,AT2> right_type;
+ typedef CML_DEFAULT_BASIS_ORIENTATION basis_orient;
+
+ /* Get matrix size: */
+ enum {
+ array_rows = left_type::array_size,
+ array_cols = right_type::array_size
+ };
+
+ /* Deduce the corresponding matrix types for the vectors: */
+ typedef CML_DEFAULT_ARRAY_LAYOUT layout;
+ typedef typename select_if<
+ array_rows == -1, dynamic<>, fixed<array_rows,1>
+ >::result left_storage;
+ typedef cml::matrix<E1,left_storage,basis_orient,layout> left_matrix;
+
+ typedef typename select_if<
+ array_cols == -1, dynamic<>, fixed<1,array_cols>
+ >::result right_storage;
+ typedef cml::matrix<E2,right_storage,basis_orient,layout> right_matrix;
+
+ /* Finally, promote the matrix types to get the result: */
+ typedef typename et::MatrixPromote<left_matrix,right_matrix>::type type;
+ typedef typename type::temporary_type temporary_type;
+ };
+
+ /** Remove const and & from the to-be-promoted types. */
+ typedef typename remove_const<
+ typename remove_reference<LeftT>::type>::type LeftBaseT;
+ typedef typename remove_const<
+ typename remove_reference<RightT>::type>::type RightBaseT;
- /* The deduced temporary type: */
+ typedef typename MatrixPromoteHelper<LeftBaseT,RightBaseT>::type type;
typedef typename type::temporary_type temporary_type;
};
{
typedef typename MatrixPromote<
typename Mat1_T::temporary_type, typename Mat2_T::temporary_type
- >::temporary_type temporary_type;
+ >::temporary_type temporary_type;
typedef typename temporary_type::value_type value_type;
};
typedef typename MatrixPromote<
typename Mat1_T::temporary_type,
typename MatrixPromote<
- typename Mat2_T::temporary_type, typename Mat3_T::temporary_type
+ typename Mat2_T::temporary_type,
+ typename Mat3_T::temporary_type
>::temporary_type
- >::temporary_type temporary_type;
+ >::temporary_type temporary_type;
typedef typename temporary_type::value_type value_type;
};
typename Mat3_T::temporary_type,
typename Mat4_T::temporary_type
>::temporary_type
- >::temporary_type
- >::temporary_type temporary_type;
+ >::temporary_type
+ >::temporary_type temporary_type;
typedef typename temporary_type::value_type value_type;
};
-/** Type promotion for a matrix and a scalar. */
-template<typename E, class AT, typename BO, typename L, typename S>
-struct MatrixPromote<cml::matrix<E,AT,BO,L>, S>
-{
- /* The deduced matrix result type (the array type is the same): */
- typedef cml::matrix<typename ScalarPromote<E,S>::type, AT, BO, L> type;
-
- /* The deduced temporary type: */
- typedef typename type::temporary_type temporary_type;
-};
-
-/** Type promotion for a scalar and a matrix. */
-template<typename S, typename E, class AT, typename BO, typename L>
-struct MatrixPromote<S, cml::matrix<E,AT,BO,L> >
-{
- /* The deduced matrix result type (the array type is the same): */
- typedef cml::matrix<typename ScalarPromote<S,E>::type, AT, BO, L> type;
-
- /* The deduced temporary type: */
- typedef typename type::temporary_type temporary_type;
-};
-
-/** Type promotion for outer product. */
-template<typename E1, class AT1, typename E2, class AT2>
-struct MatrixPromote< cml::vector<E1,AT1>, cml::vector<E2,AT2> >
-{
- typedef cml::vector<E1,AT1> left_type;
- typedef cml::vector<E2,AT2> right_type;
- typedef CML_DEFAULT_BASIS_ORIENTATION basis_orient;
-
- /* Get matrix size: */
- enum {
- array_rows = left_type::array_size,
- array_cols = right_type::array_size
- };
-
- /* Deduce the corresponding matrix types for the vectors: */
- typedef CML_DEFAULT_ARRAY_LAYOUT layout;
- typedef typename select_if<
- array_rows == -1, dynamic<>, fixed<array_rows,1>
- >::result left_storage;
- typedef cml::matrix<E1,left_storage,basis_orient,layout> left_matrix;
-
- typedef typename select_if<
- array_cols == -1, dynamic<>, fixed<1,array_cols>
- >::result right_storage;
- typedef cml::matrix<E2,right_storage,basis_orient,layout> right_matrix;
-
- /* Finally, promote the matrix types to get the result: */
- typedef typename et::MatrixPromote<left_matrix,right_matrix>::type type;
- typedef typename type::temporary_type temporary_type;
-};
-
} // namespace et
} // namespace cml
private:
/* XXX Blah, a temp. hack to fix the auto-resizing stuff below. */
- matrix_size hack_actual_size(const SrcT& /*src*/, scalar_result_tag) {
- return matrix_size(1,1);
+ matrix_size hack_actual_size(
+ matrix_type& dest, const SrcT& /*src*/, scalar_result_tag
+ )
+ {
+ typedef ExprTraits<matrix_type> dest_traits;
+ return dest_traits().size(dest);
}
- matrix_size hack_actual_size(const SrcT& src, matrix_result_tag) {
+ matrix_size hack_actual_size(
+ matrix_type& /*dest*/, const SrcT& src, matrix_result_tag
+ )
+ {
typedef ExprTraits<SrcT> src_traits;
return src_traits().size(src);
}
#if defined(CML_AUTOMATIC_MATRIX_RESIZE_ON_ASSIGNMENT)
/* Get the size of src. This also causes src to check its size: */
matrix_size N = hack_actual_size(
- src, typename src_traits::result_tag());
+ dest, src, typename src_traits::result_tag());
/* Set the destination matrix's size: */
dest.resize(N.first,N.second);
{
return CheckedSize(dest,src,dynamic_size_tag());
}
- /* XXX Blah, a temp. hack to fix the auto-resizing stuff below. */
+
+
public:
return (*this = *this * e);
}
+ /** Return access to the data as a raw pointer. */
+ typename vector_type::pointer data() { return m_q.data(); }
+
+ /** Return access to the data as a const raw pointer. */
+ const typename vector_type::pointer data() const { return m_q.data(); }
+
+
/* NOTE: Quaternion division no longer supported, but I'm leaving the
code here for reference (Jesse) */
#include <cstdlib> // For std::rand.
#include <cml/constants.h>
+#if defined(_MSC_VER)
+#pragma push_macro("min")
+#pragma push_macro("max")
+#undef min
+#undef max
+#endif
+
namespace cml {
/** Sign of input value as double. */
/** Test input value for inclusion in [min, max]. */
template < typename T >
bool in_range(T value, T min, T max) {
- return value >= min && value <= max;
+ return !(value < min) && !(value > max);
}
/** Map input value from [min1, max1] to [min2, max2]. */
}
-/** For convenient squaring of expressions. */
+/** Square a value. */
template < typename T >
T sqr(T value) {
return value * value;
}
+/** Cube a value. */
+template < typename T >
+T cub(T value) {
+ return value * value * value;
+}
+
/** Inverse square root. */
template < typename T >
T inv_sqrt(T value) {
return std::sqrt(length_squared(x,y,z));
}
+/** Index of maximum of 2 values. */
+template < typename T >
+size_t index_of_max(T a, T b) {
+ return a > b ? 0 : 1;
+}
+
+/** Index of maximum of 2 values by magnitude. */
+template < typename T >
+size_t index_of_max_abs(T a, T b) {
+ return index_of_max(std::fabs(a),std::fabs(b));
+}
+
+/** Index of minimum of 2 values. */
+template < typename T >
+size_t index_of_min(T a, T b) {
+ return a < b ? 0 : 1;
+}
+
+/** Index of minimum of 2 values by magnitude. */
+template < typename T >
+size_t index_of_min_abs(T a, T b) {
+ return index_of_min(std::fabs(a),std::fabs(b));
+}
+
/** Index of maximum of 3 values. */
template < typename T >
size_t index_of_max(T a, T b, T c) {
} // namespace cml
+#if defined(_MSC_VER)
+#pragma pop_macro("min")
+#pragma pop_macro("max")
+#endif
+
#endif
// -------------------------------------------------------------------------
/* Shorthand for the type of this vector: */
typedef vector<Element,generator_type> vector_type;
+ /* The vector coordinate type: */
+ typedef Element coordinate_type;
+
/* For integration into the expression template code: */
typedef vector_type expr_type;
/* Shorthand for the type of this vector: */
typedef vector<Element,generator_type> vector_type;
+ /* The vector coordinate type: */
+ typedef Element coordinate_type;
+
/* For integration into the expression template code: */
typedef vector_type expr_type;
typedef cml::et::assignable_tag assignable_tag;
+ public:
+
+ /** Static constant containing the vector's space dimension. */
+ enum { dimension = Size };
+
+
public:
/** Return square of the length. */
/* Shorthand for the type of this vector: */
typedef vector<Element,generator_type> vector_type;
+ /* The vector coordinate type: */
+ typedef Element coordinate_type;
+
/* For integration into the expression template code: */
typedef vector_type expr_type;
typedef cml::et::assignable_tag assignable_tag;
+ public:
+
+ /** Static constant containing the vector's space dimension. */
+ enum { dimension = Size };
+
+
public:
/** Return square of the length. */
}
}
+ /** Return a subvector by removing element i.
+ *
+ * @internal This is horribly inefficient...
+ */
+ subvector_type subvector(size_t i) const {
+ subvector_type s;
+ for(size_t m = 0, n = 0; m < this->size(); ++ m)
+ if(m != i) s[n++] = (*this)[m];
+ return s;
+ };
+
public:
template<typename E, class AT > inline std::ostream&
operator<<(std::ostream& os, const vector<E,AT>& v)
{
- os << "[";
for (size_t i = 0; i < v.size(); ++i) {
os << " " << v[i];
}
- os << " ]";
return os;
}
template< class XprT > inline std::ostream&
operator<<(std::ostream& os, const et::VectorXpr<XprT>& v)
{
- os << "[";
for (size_t i = 0; i < v.size(); ++i) {
os << " " << v[i];
}
- os << " ]";
return os;
}
projects / chaz / yoink / commitdiff
