What about my vector implementation?

Started by
19 comments, last by cignox1 16 years, 6 months ago
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
Advertisement
Quote:
//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); }

I find this to be incredibly annoying, especially your overload of operator|

Named functions make the operations signifigantly more understandable.
Quote:Original post by jpetrie
Quote:
//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); }

I find this to be incredibly annoying, especially your overload of operator|

Named functions make the operations signifigantly more understandable.


Yes, in my old implementation I had both the operators overloading and the named functions. I thought to switch to them, not yet decided though. I will definitely add the methods, anyway... Anything else? Thak you for answering!

EDIT: I will clarify my embarassingly use of the | operator: I've always used * and ^ (* looks like a dot, ^ remembers the right angle between two orthogonal vectors). I saw someone using * and | over the net, so when I thought to have a 3d specialized dot product, I took that operator. Not my best hour, I suppose :-)
For dot and cross products I tend to do something like:
template <typename T, typename Op>struct LtProxy {  public:    LtProxy(T r) : ref(r) {}    T ref;  private:    LtProxy & operator=(const LtProxy &);};const struct dot_t {} dot = {};const struct cross_t {} cross = {};inline LtProxy<const Vector &, dot_t> operator<(const Vector & lhs, dot_t) {  return LtProxy<const Vector &, dot_t>(lhs);}inline LtProxy<const Vector &, cross_t> operator<(const Vector & lhs, cross_t) {  return LtProxy<const Vector &, cross_t>(lhs);}inline double operator>(const LtProxy<const Vector &, dot_t> & lhs, const Vector & rhs) {  return dot_product(lhs.ref, rhs);}inline Vector operator>(const LtProxy<const Vector &, cross_t> & lhs, const Vector & rhs) {  return cross_product(lhs.ref, rhs);}

So then I can use "A <dot> B" or "A <cross> B" to calculate the dot or cross product. On the down side, it's not idiomatic C++, so it tends to generate WTF comments from people seeing it used before seeing how its implemented.
Quote:Original post by SiCrane
For dot and cross products I tend to do something like:
*** Source Snippet Removed ***
So then I can use "A <dot> B" or "A <cross> B" to calculate the dot or cross product. On the down side, it's not idiomatic C++, so it tends to generate WTF comments from people seeing it used before seeing how its implemented.


WTF?
just a joke :-) Interesting use of templates, reminds to me some uses of Ocaml (well, a sort of...)
I don't think I will use them, but thank you anyway. By the way, nobody made comments about the overall use of my template... should I assume that it is a acceptable vector4 class (personal tastes aside?) :-)
The use of templates is perfectly acceptable and is the way I implemented a vector class. The only thing I did different was the naming:
template <T>class Vector4Base{  // stuff};typedef Vector4Base<float> Vector4;

which cut down on the typing when using the vector type and allowed for a global change to the instantiated type with one line of source changed.

Redefining the meaning of operators is generally frowned upon (i.e. using * for a cross product) as it requires knowing the implementation to understand the code, a third party looking at the code will assume the standard definition for the operator, thus misinterpreting the code. Always use a named method for undefined operations (i.e. not *,/,+,-,etc).

I like SiCrane's solution. But then again I do like to write code that confuses the hell out of people:
(expr ? fn1 : fn2) (args);

was a gem I came up with many years ago (it's a code size optimisation, honest). Duff's device is truely awesome. I'm currently abusing .Net's reflection system to create GUI layout editors.

Skizz
Quote:Original post by SiCrane
For dot and cross products I tend to do something like:
*** Source Snippet Removed ***
So then I can use "A <dot> B" or "A <cross> B" to calculate the dot or cross product. On the down side, it's not idiomatic C++, so it tends to generate WTF comments from people seeing it used before seeing how its implemented.

And I thought my math library was over-engineered [wink]
Quote:Original post by Skizz
The use of templates is perfectly acceptable and is the way I implemented a vector class.


Thank you, actually I was meaning my implementation (major issues I'm understimating, wrong design approaches and so on) :-)

Quote:
The only thing I did different was the naming:
template <T>class Vector4Base{  // stuff};typedef Vector4Base<float> Vector4;

which cut down on the typing when using the vector type and allowed for a global change to the instantiated type with one line of source changed.


Don't worry, I will typedef the class for the most common types (floats, doubles and ints).

Quote:
Redefining the meaning of operators is generally frowned upon (i.e. using * for a cross product) as it requires knowing the implementation to understand the code, a third party looking at the code will assume the standard definition for the operator, thus misinterpreting the code. Always use a named method for undefined operations (i.e. not *,/,+,-,etc).


Well, apparently nobody like overloaded operators (at least used in this manner). Most probably I will move to named functions then.
I like to make operator *(const Vector &, const Vector &) an element-wise multiplication, which is useful for things like colour combining. It also means that vector * 4 will have the same result as vector * Vector(4).

I mucked around with templating the hell out of my vector classes at one point, but in the end I found that keeping Vector3f, Vector3d and Vector3i separate gave me far more flexibility in the implementation, and debugging was far easier. The "A <dot> B" and "A <cross> B" thing is clever, but again it just makes things more confusing for anyone else trying to understand your code and complicates debugging. I prescribe to a.dot(b), but some people like dot(a, b) because it's more like a (operand neutral) maths function.
Quote:Original post by SiCrane
For dot and cross products I tend to do something like:
*** Source Snippet Removed ***
So then I can use "A <dot> B" or "A <cross> B" to calculate the dot or cross product. On the down side, it's not idiomatic C++, so it tends to generate WTF comments from people seeing it used before seeing how its implemented.


That is one of the most ridiculous things I have ever seen. It deserves to be submitted to The Daily WTF. I am really surprised that someone of your caliber would write that. Oh well, I guess everyone has their vices. It should be obvious that I am not a fan of Boost.Spirit.
John BoltonLocomotive Games (THQ)Current Project: Destroy All Humans (Wii). IN STORES NOW!

This topic is closed to new replies.

Advertisement