# What about my vector implementation?

## Recommended Posts

cignox1    735
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
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) {}

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]); }

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);*/

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 on other sites
jpetrie    13103
Quote:
 //Dot, Dot3 and Cross productsconst T operator * (const vector4& v) const { return (x*x + y*y + z*z + w*w); }const T operator | (const vector4& v) const { return (x*x + y*y + z*z); }const vector4 operator ^ (const vector4& v) const { return vector4(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 on other sites
cignox1    735
Quote:
Original post by jpetrie
Quote:
 //Dot, Dot3 and Cross productsconst T operator * (const vector4& v) const { return (x*x + y*y + z*z + w*w); }const T operator | (const vector4& v) const { return (x*x + y*y + z*z); }const vector4 operator ^ (const vector4& v) const { return vector4(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 on other sites
SiCrane    11839
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 on other sites
cignox1    735
Quote:
 Original post by SiCraneFor dot and cross products I tend to do something like:*** Source Snippet Removed ***So then I can use "A B" or "A 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 on other sites
Skizz    794
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 on other sites
Zipster    2359
Quote:
 Original post by SiCraneFor dot and cross products I tend to do something like:*** Source Snippet Removed ***So then I can use "A B" or "A 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 on other sites
cignox1    735
Quote:
 Original post by SkizzThe 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 class Vector4Base{ // stuff};typedef Vector4Base 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 on other sites
hh10k    589
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 on other sites
JohnBolton    1372
Quote:
 Original post by SiCraneFor dot and cross products I tend to do something like:*** Source Snippet Removed ***So then I can use "A B" or "A 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 on other sites
Hodgman    51223
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 SiCraneFor dot and cross products I tend to do something like:*** Source Snippet Removed ***So then I can use "A B" or "A 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 on other sites
cignox1    735
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 on other sites
ToohrVyk    1595
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 on other sites
rozz666    896
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 on other sites
cignox1    735
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 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 on other sites
rozz666    896
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 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 on other sites
SiCrane    11839
Quote:
 Original post by JohnBoltonThat 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 on other sites
d00fus    328
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 on other sites
cignox1    735
Quote:
 Original post by d00fusToohrVyk 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 on other sites
clb    2147

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 on other sites
cignox1    735
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!