Sign in to follow this  

What about my vector implementation?

This topic is 3716 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

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[i]); }
		const T& operator[](std::size_t i) const { return (*this).*(vec_elem_ptrs[i]); }

		//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

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites
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 :-)

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites
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?) :-)

Share this post


Link to post
Share on other sites
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

Share this post


Link to post
Share on other sites
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]

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites
Yeah, I like to have * just do element-wise multiplication too, and have named functions for dot and cross multiplication (because everyone seems to have a differing opinion on which operators to use). Only use operator overloading if it's going to be obvious what the operator does - If I saw vec1 = vec2|vec3, I would assume it was some kind of element-wise logical-or...

Quote:
Original post by JohnBolton
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.

Ridiculous? Yes.
Confusing? Yes.
Smothered in awesome sauce? YES!!

Share this post


Link to post
Share on other sites
I see many responses... cool :-)
But I must say I'm a bit surprised about the fact that almost all the replies are about the use of operator overloading for dot and cross products. This is a topic covered in many places and I'm aware of the fact that operators overloading can lead to confusing practices and should be used only where overloaded operators have a behaviour semantically similar to the original operators (that's why my previous class also had named methods to perform dot and cross products: I just did't add them to the new version yet).

However, I'm more interested in your advices about the class as a whole, i.e how I used templates. const correctness, mathematical operations, interface design and so on (I suppose that nothing really wrong was made, otherwise you would have pointed it out :-)...
Com'on guys, Do you have some more observations not about those operators? ;-)

QUESTION:
When you need to add a scalar value (i.e. 0 or 1) which form do you use?
-T(0)
-(T)0
-T(0.0)
-(T)0.0

And are there reasons to choose one over the others?

This doesn't mean I'm not grateful for your replies, of course!

Thank you again!
Alessandro

Share this post


Link to post
Share on other sites
You need a typedef T element_type; in your class to make generic programmer easier (so that you can deduce what T is in a template function where the vector is a parameter).

Share this post


Link to post
Share on other sites
cignox1, you may consider using expression templates so that code like:
Vector4<float> v5 = v1 + v2 + v3 + v4;

doesn't produce unneccesssary temporary objects.

Share this post


Link to post
Share on other sites
Quote:
You need a typedef T element_type; in your class to make generic programmer easier (so that you can deduce what T is in a template function where the vector is a parameter).


Quote:

cignox1, you may consider using expression templates so that code like:
Vector4<float> v5 = v1 + v2 + v3 + v4;

doesn't produce unneccesssary temporary objects.


Could you both explain to me what you mean? I'm not sure to fully understand...

EDIT: perhaps I understood what ToohrVyk meant: I start to implement the ray class. I will parametrize it with the vector type (vector4 or vector2 will be accepted):
template<typename T>
class ray
{
public:
T origin, direction;
};

If I put the typedef inside the vector class I can then do something like this:
T.element_type a = ...
where I need temp values for calculations, with the correct type without requiring to ask for a second type parameter when instantiating the template. Did I get it?

[Edited by - cignox1 on October 12, 2007 7:55:20 AM]

Share this post


Link to post
Share on other sites
Quote:
Original post by cignox1
Quote:
You need a typedef T element_type; in your class to make generic programmer easier (so that you can deduce what T is in a template function where the vector is a parameter).


Quote:

cignox1, you may consider using expression templates so that code like:
Vector4<float> v5 = v1 + v2 + v3 + v4;

doesn't produce unneccesssary temporary objects.


Could you both explain to me what you mean? I'm not sure to fully understand...


T is a typename, but this does notimply that this is a type of vector elements.
Making a typedef element_type will make it clear that it is.

Expression templates are a great way to speed up your code. You can read about it here or just google gamedev for it.

Share this post


Link to post
Share on other sites
Quote:
Original post by JohnBolton
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.


If it's any consolation, the header that implements that in my code library is named "evil.h". It's not really Daily WTF worthy though, since I've never used it in anything but personal projects.

Share this post


Link to post
Share on other sites
ToohrVyk meant that in the public section of your vector class, add


template<typename T>
class vector4
{
public:
typedef T element_type; // add this
...
};



This means you can deduce the element type from a vector passed as a parameter to a function template, as in


template <typename T>
void DoSomethingWithVector(const T& vec)
{
T::element_type val = vec.Foo();
...
}



You can implement generic algorithms in this way that aren't dependent on the type vector4 was instantiated with.

Share this post


Link to post
Share on other sites
Quote:
Original post by d00fus
ToohrVyk meant that in the public section of your vector class, add

*** Source Snippet Removed ***

This means you can deduce the element type from a vector passed as a parameter to a function template, as in

*** Source Snippet Removed ***

You can implement generic algorithms in this way that aren't dependent on the type vector4 was instantiated with.



Yep, I realized that just after asking for explanation :-) In fact I added a few words in my previous post to get confirm to my supposition. Incidentally, you reminded me one more time that access to static field in c++ is NOT the same as in c# or Java ;-)

This is a good idea, one of those things that make me wonder if C++ has an actually finite amount of features. After years of programming, one still sees new things. SiCrane post is another one, though I suppose that templates add a new whole set of features...

Share this post


Link to post
Share on other sites
(I'm not sure if you already know about this, so ignore me if I sound patronizing)

Before you bash on to implementing your vector4 for use in a raytracer, I'd like to note that having an SSE-enabled vector4 structure is not that useful. To get most out of SIMD, you'll probably want to have a SSE-enabled structure that stores 4 3D vectors (12 elements) in a SoA (structure of arrays) manner. For example, see Using Intel Streaming SIMD Extensions for 3D Geometry Processing.

Share this post


Link to post
Share on other sites
Quote:
Original post by clb
(I'm not sure if you already know about this, so ignore me if I sound patronizing)

Before you bash on to implementing your vector4 for use in a raytracer, I'd like to note that having an SSE-enabled vector4 structure is not that useful. To get most out of SIMD, you'll probably want to have a SSE-enabled structure that stores 4 3D vectors (12 elements) in a SoA (structure of arrays) manner. For example, see Using Intel Streaming SIMD Extensions for 3D Geometry Processing.


Yes, SSE can for sure be used better than just for matrix multiplication. It could also help ray packing AFAIK. I'm not going to make realtime RT though, so the architecture I want to build is the standard (and easy) one. In addition, I've never done much low level optimization, and I'm sure that what I would gain from using SSE better would be wasted elsewere, so I'm going to only speed up those math routines that don't require design changes :-)
Thank you anyway!

Share this post


Link to post
Share on other sites

This topic is 3716 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this