Hi, all! I've posted a couple of weeks ago to have tips about how to implement my new geometry classes. I've moved to templates and I will specialize it for SSE. I would like to know what do you think about my first class (vector4): are templates used in the correct way? Did I skip some important operation?
Thank to those willing to share their experience.
Note that this class will be used for my next raytracer, so it's not going to be used for generic mathematics problems.
#ifndef GEOMETRY_H
#define GEOMETRY_H
#include <cmath>
#include <iostream>
namespace qlgeom
{
template<typename T> class vector4;
template<typename T> class point4;
template<typename T> class vector2;
template<typename T> class point2;
template<typename T> class matrix4x4;
template<typename T> class matrix2x2;
template<typename T> class ray;
template<typename T> class quaternion;
/*
The vector4 class defines a templated 4-dimensional vector. Note that in order to use it for homogeneous coordinates
the user must keep attention not to mix points and vectors (that is w = 1 and w = 0) where they are not supposed to be
used. For example, Point - Point works (and results in a vector), while Vector - Vector does not. Length() gives the
Length() of the 4 dimensions vector, not the 3d vector in the homogeneous coordinate space. This means that taking the
length of a point (w = 1) will not produce the desidered result in most cases. The same for most of the other operators
and methods.
Additional utility methods are provided to help in the use with homogeneous coordinate systems. These methods will
correctly address the w coordinate.
Remember that vectors and points have different semantics and that even while both them can be expressed as 4d vectors
of elements, standard operations that doesn't address those issues expliticy will work with homogeneous coordinates only
if the user takes care about using vectors and points in the right places and about checking the value of w.
*/
template<typename T>
class vector4
{
public:
T x, y, z, w;
vector4() {}
vector4(const T s) : x(s), y(s), z(s), w(s) {}
vector4(const T _x, const T _y, const T _z) : x(_x), y(_y), z(_z), w(0) {}
vector4(const T _x, const T _y, const T _z, const T _w) : x(_x), y(_y), z(_z), w(_w){}
vector4(const vector4<T>& src) : x(src.x), y(src.y), z(src.z), w(src.w) {}
//Add indexed access capabilities
T& operator[](std::size_t i) { return (*this).*(vec_elem_ptrs); }
const T& operator[](std::size_t i) const { return (*this).*(vec_elem_ptrs); }
//Overloaded operators
friend const vector4<T> operator + (const T l, const vector4<T>& r);
friend const vector4<T> operator - (const T l, const vector4<T>& r);
friend const vector4<T> operator * (const T l, const vector4<T>& r);
friend const vector4<T> operator / (const T l, const vector4<T>& r);
bool operator == (const vector4<T>& v) const { return (x == v.x) && (y == v.y) && (z == v.z) && (w == v.w); }
bool operator != (const vector4<T>& v) const { return !(*this == v); }
const vector4<T>& operator = (const vector4<T>& v) { x = v.x; y = v.y; z = v.z; w = v.w; return *this; }
const vector4<T> operator + (const vector4<T>& v) const { return vector4<T>(x+v.x, y+v.y, z+v.z, w+v.w); }
const vector4<T> operator - (const vector4<T>& v) const { return vector4<T>(x-v.x, y-v.y, z-v.z, w-v.w); }
const vector4<T>& operator += (const vector4<T>& v) { x += v.x; y += v.y; z += v.z; w += v.w; return *this; }
const vector4<T>& operator -= (const vector4<T>& v) { x -= v.x; y -= v.y; z -= v.z; w -= v.w; return *this; }
const vector4<T> operator * (const T s) const { return vector4<T>(x*s, y*s, z*s, w*s); }
const vector4<T> operator / (const T s) const { return vector4<T>(x/s, y/s, z/s, w/s); }
const vector4<T>& operator *= (const T s) const { x *= s; y *= s; z *= s; w *= s; return *this; }
const vector4<T>& operator /= (const T s) const { x /= s; y /= s; z /= s; w /= s; return *this; }
//Dot, Dot3 and Cross products
const T operator * (const vector4<T>& v) const { return (x*x + y*y + z*z + w*w); }
const T operator | (const vector4<T>& v) const { return (x*x + y*y + z*z); }
const vector4<T> operator ^ (const vector4<T>& v) const { return vector4<T>(y*v.z - z*v.y, x*v.z - z*v.x, x*v.y - y*v.z, (T)0); }
//Utility methods
void Normalize()
{
const T t = T(1) / Length();
x *= t; y *= t; z *= t; w *= t;
}
const T Length() const { return T(Sqrt(x*x + y*y + z*z + w*w)); }
const T InvLength() const { return T(1) / Length(); }
const T Magnitude() const { return (T) (x*x + y*y + z*z + w*w); }
int DominantAxisInd3() const { return (x > y) ? ((x > z) ? 0 : 2) : ((y > z) ? 1 : 2); }
const T DominantAxis3() const { return (x > y) ? ((x > z) ? x : z) : ((y > z) ? y : z); }
//Checking for non zero vector is left to the caller!
const vector4<T> Direction(void) const
{
vector4<T> t(*this);
t.Normalize();
return t;
}
void Invert()
{
x = -x; y = -y; z = -z; w = -w;
}
const vector4<T> Project(const vector4<T>& v) const
{
const T ll = Length();
return ((*this * src) / (ll * ll)) * v;
}
T Angle(const vector4<T> &v) const
{
return Acos(*this * v)/(Length() * v.Length());
}
const vector4<T> Reflect(const vector4<T> &n) const
{
return *this - (T(2.0) * ( *this * n)) * n;
}
const vector4 Refract(const vector4<T>& n, T ior_i, T ior_r, T& tir) const
{
// calculate refraction
const T rior = ior_i / ior_r;
const T cos_i = -(n * *this);
const T cos_t = T(1) - rior * rior * (T(1) - cos_i * cos_i);
if (cos_t > T(0))
{
vector4 res((rior * *this) + (rior * cos_i - Sqrt( cos_t )) * n);
tir = false;
return res;
}
else tir = true;
return Reflect(n);
}
//Calculates the linear interpolation in the segment given by [this, end] at the distance x = (0,1)
void Lerp(const vector4<T>& v1, const vector4<T>& v2, const float t)
{
*this = v1 + (v2 - v1) * t;
}
//Calculates the linear interpolation of the angle between this and v in the point x = (0,1)
void Slerp(const vector4<T>& v1, const vector4<T>& v2, const float t)
{
SlerpUnit(v1.Direction(), v2.Direction(), t);
}
//Calcuales the linear interpolation of the angle between this and v in the point x = (0,1)
//Both *this and v must be unit length!
void SlerpUnit(const vector4<T>& v1, const vector4<T>& v2, const float t)
{
if ( t <= T(0.0) )
{
(*this) = v1;
return;
}
else if ( t >= T(1.0) )
{
(*this) = v2;
return;
}
const T cosom = v1 * v2;
T scale0, scale1;
if ( ( T(1.0) - cosom ) > Q_MATCH_FACTOR )
{
const T omega = Acos( cosom );
const T sinom = Sin( omega );
scale0 = Sin( ( T(1.0) - t ) * omega ) / sinom;
scale1 = Sin( t * omega ) / sinom;
}
else
{
scale0 = T(1.0) - t;
scale1 = t;
}
(*this) = ( v1 * scale0 + v2 * scale1 );
Normalize();
}
//Help methods for Homogeneous coordinates
void Reduce(void)
{
if(w == T(0) || w == T(1)) return;
x /= w;
y /= w;
z /= w;
w = T(1.0);
}
static void CoordinateSystem3(const vector4<T>& v1, vector4<T>& v2, vector4<T>& v3);
T Volume3(const vector4<T>& v1, const vector4<T>& v2)
{
return Abs(((*this ^ v1) * v2));
}
T Area3(const vector4<T>& v)
{
return (*this ^ v).Length();
}
//Debug methods
private:
typedef T vector4<T>::* const vec_elem_ptrs_t[4];
static vec_elem_ptrs_t vec_elem_ptrs;
};
//Initialize the static pointer-to-member variable
template<typename T>
typename vector4<T>::vec_elem_ptrs_t vector4<T>::vec_elem_ptrs = { &vector4<T>::x, &vector4<T>::y, &vector4<T>::z, &vector4<T>::w };
};
#endif //GEOMETRY_H
#include "./geometry.h"
namespace qlgeom
{
/*----------------Implementation of the vector4 class---------------------------------------*/
template <typename T>
void vector4<T>::CoordinateSystem3(const vector4<T>& v1, vector4<T>& v2, vector4<T>& v3)
{
if(Abs(v1.x) > Abs(v1.y))
{
T invLen = T(1) / Sqrt(v1.x * v1.x + v1.z * v1.z);
v2 = vector4(-v1.z * invLen, T(0), v1.x * invLen);
}
else
{
T invLen = T(1.0) / Sqrt(v1.y * v1.y + v1.z * v1.z);
v2 = vector4(T(0), v1.z * invLen, -v1.y * invLen);
}
v3 = v1 ^ v2;
}
/* friend const vector4<T> operator + (const T l, const vector4<T>& r);
friend const vector4<T> operator - (const T l, const vector4<T>& r);
friend const vector4<T> operator * (const T l, const vector4<T>& r);
friend const vector4<T> operator / (const T l, const vector4<T>& r);*/
//More operator overloading
template <typename T>
const vector4<T> operator + (T l, const vector4<T>& r)
{
return vector4<T>(r.x + l, r.y + l, r.z + l);
}
template <typename T>
const vector4<T> operator - (T l, const vector4<T>& r)
{
return vector4<T>(r.x - l, r.y - l, r.z - l);
}
template <typename T>
const vector4<T> operator * (T l, const vector4<T>& r)
{
return vector4<T>(r.x * l, r.y * l, r.z * l);
}
template <typename T>
const vector4<T> operator / (T l, const vector4<T>& r)
{
T div = T(1.0) / l; return vector4<T>(r.x * div, r.y * div, r.z * div);
}
}
My next step will be to rewrite my matrix and ray classes using templates.
Let me know!
Alessandro
EDIT: the 'real' type was not even defined (was typedefed in my previous class), don't know why the compiler did't complain :-( I changed it to T