Sign in to follow this  

Is my matrix math busted - it probably is?

This topic is 416 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

I am fighting with my 2x2 and 4x4 rotation matrices right now, havent written any math in years...
 
When i define a normal from an axis with 90 degrees i get a proper value of cos = 0 and sin = 1 which is correct - pointing upwards in a cartesian coordinate system, moving counter clock wise when the radians increase.
However in my physics system am working on, i am getting the opposite: Clockwise rotation... so a plane defined with half a pi points downwards O_o
 
Are my math totally busted or are there just some simple mistake i made - i dont see it.
 
typedef float F32;
typedef uint32_t U32;

inline F32 Sin(F32 radians) {
	F32 result = sinf(radians);
	return(result);
}
inline F32 Cos(F32 radians) {
	F32 result = cosf(radians);
	return(result);
}
inline F32 ATan2(F32 y, F32 x)
{
	F32 result = atan2f(y, x);
	return(result);
}

union Vec2f {
	struct {
		F32 x, y;
	};
	struct {
		F32 u, v;
	};
	F32 p[2];
};
union Vec4f {
	struct {
		union {
			Vec3f xyz;
			struct {
				F32 x, y, z;
			};
		};
		F32 w;
	};
	struct {
		union {
			Vec3f rgb;
			struct {
				F32 r, g, b;
			};
		};
		F32 a;
	};
	struct {
		Vec3f xyz;
		F32 w;
	};
	struct
	{
		Vec2f xy;
		F32 ignored0;
		F32 ignored1;
	};
	struct
	{
		F32 ignored2;
		Vec2f yz;
		F32 ignored3;
	};
	struct
	{
		F32 ignored4;
		F32 ignored5;
		Vec2f zw;
	};
	F32 p[4];
};
union Mat2f {
	struct {
		Vec2f col1;
		Vec2f col2;
	};
	F32 m[4];
};
union Mat4f {
	struct {
		Vec4f col1;
		Vec4f col2;
		Vec4f col3;
		Vec4f col4;
	};
	F32 m[16];
};

inline Vec2f V2(F32 x, F32 y) {
	Vec2f result;
	result.x = x;
	result.y = y;
	return(result);
}

inline F32 Dot(const Vec2f &a, const Vec2f &b)
{
	F32 result = a.x * b.x + a.y * b.y;
	return(result);
}

inline Vec2f operator *(const Vec2f &a, const Mat2f &mat)
{
	Vec2f result = V2(Dot(a, mat.col1), Dot(a, mat.col2));
	return(result);
}

inline Mat2f MakeIdentity2() {
	Mat2f result = {};
	result.col1 = V2(1.0f, 0.0f);
	result.col2 = V2(0.0f, 1.0f);
	return (result);
}

inline Mat2f MakeRotation2(F32 angle) {
	F32 s = Sin(angle);
	F32 c = Cos(angle);
	Mat2f result = {};
	result.col1 = V2(c, s);
	result.col2 = V2(-s, c);
	return(result);
}

inline Mat2f Transpose2(const Mat2f &m) {
	Mat2f result = {};
	result.col1 = V2(m.col1.x, m.col2.x);
	result.col2 = V2(m.col1.y, m.col2.y);
	return(result);
}

inline Vec2f operator *(const Mat2f &a, const Vec2f &v) {
	Vec2f result = V2(Dot(a.col1, v), Dot(a.col2, v));
	return(result);
}

inline Mat2f operator *(const Mat2f &a, const Mat2f &b)
{
	Mat2f result = {};
	result.col1 = a * b.col1;
	result.col2 = a * b.col2;
	return(result);
}

inline F32 ToAngle(const Mat2f &mat) {
	F32 result = ATan2(mat.col1.y, mat.col1.x);
	return(result);
}

inline Mat4f MakeIdentity4() {
	Mat4f result = {};
	result.col1 = V4(1.0f, 0.0f, 0.0f, 0.0f);
	result.col2 = V4(0.0f, 1.0f, 0.0f, 0.0f);
	result.col3 = V4(0.0f, 0.0f, 1.0f, 0.0f);
	result.col4 = V4(0.0f, 0.0f, 0.0f, 1.0f);
	return (result);
}

inline Mat4f MakeRotation4(const Mat2f &mat2) {
	Mat4f result = MakeIdentity4();
	result.col1.xy = mat2.col1;
	result.col2.xy = mat2.col2;
	return (result);
}

inline Mat4f MakeTranslation4(const Vec2f &p) {
	Mat4f result = MakeIdentity4();
	result.col4.x = p.x;
	result.col4.y = p.y;
	result.col4.z = 0.0f;
	return (result);
}

inline Mat4f operator *(const Mat4f &a, const Mat4f &b)
{
	// http://stackoverflow.com/questions/18499971/efficient-4x4-matrix-multiplication-c-vs-assembly
	Mat4f result = MakeIdentity4();
	for (U32 i = 0; i < 16; i += 4) {
		for (U32 j = 0; j < 4; ++j) {
			result.m[i + j] =
				(b.m[i + 0] * a.m[j + 0])
				+ (b.m[i + 1] * a.m[j + 4])
				+ (b.m[i + 2] * a.m[j + 8])
				+ (b.m[i + 3] * a.m[j + 12]);
		}
	}
	return(result);
}[size=4]
[/size]
Edited by Finalspace

Share this post


Link to post
Share on other sites
Found the culprit, the 2x2 mult was busted:

This is totally wrong - what was i thinking:
inline Mat2f operator *(const Mat2f &a, const Mat2f &b)
{
	Mat2f result = {};
	result.col1 = a * b.col1;
	result.col2 = a * b.col2;
	return(result);
}

This is correct:
inline Mat2f operator *(const Mat2f &a, const Mat2f &b)
{
	Mat2f result = {};
	for (U32 i = 0; i < 2; ++i) {
		for (U32 j = 0; j < 2; ++j) {
			F32 sum = 0;
			for (U32 k = 0; k < 2; ++k) {
				sum += a.m[i * 2 + k] * b.m[k * 2 + j];
			}
			result.m[i * 2 + j] = sum;
		}
	}
	return(result);
}
Edited by Finalspace

Share this post


Link to post
Share on other sites

Found the culprit, the 2x2 mult was busted:

This is totally wrong - what was i thinking:

inline Mat2f operator *(const Mat2f &a, const Mat2f &b)
{
	Mat2f result = {};
	result.col1 = a * b.col1;
	result.col2 = a * b.col2;
	return(result);
}

This is correct:
inline Mat2f operator *(const Mat2f &a, const Mat2f &b)
{
	Mat2f result = {};
	for (U32 i = 0; i < 2; ++i) {
		for (U32 j = 0; j < 2; ++j) {
			F32 sum = 0;
			for (U32 k = 0; k < 2; ++k) {
				sum += a.m[i * 2 + k] * b.m[k * 2 + j];
			}
			result.m[i * 2 + j] = sum;
		}
	}
	return(result);
}

 

You might want to check your Matrix*vector code, your original code was correct. By rewriting it you just bypassed a bug in your Matrix*vector code. Your M*V code should dot the rows of the matrix with the vector, not the columns. My guess is that you just copied the V*M code (which is the same as transpose(M)*V), without swapping rows/columns.

Edited by Aressera

Share this post


Link to post
Share on other sites

 

Found the culprit, the 2x2 mult was busted:

This is totally wrong - what was i thinking:

inline Mat2f operator *(const Mat2f &a, const Mat2f &b)
{
	Mat2f result = {};
	result.col1 = a * b.col1;
	result.col2 = a * b.col2;
	return(result);
}

This is correct:
inline Mat2f operator *(const Mat2f &a, const Mat2f &b)
{
	Mat2f result = {};
	for (U32 i = 0; i < 2; ++i) {
		for (U32 j = 0; j < 2; ++j) {
			F32 sum = 0;
			for (U32 k = 0; k < 2; ++k) {
				sum += a.m[i * 2 + k] * b.m[k * 2 + j];
			}
			result.m[i * 2 + j] = sum;
		}
	}
	return(result);
}

 

You might want to check your Matrix*vector code, your original code was correct. By rewriting it you just bypassed a bug in your Matrix*vector code. Your M*V code should dot the rows of the matrix with the vector, not the columns. My guess is that you just copied the V*M code (which is the same as transpose(M)*V), without swapping rows/columns.

 

 

Oh thanks, you are right. The matrix * vector was wrong and the multiplication was correct. Fixed it. Need to add a unit test for this - so this will never happen again.

Share this post


Link to post
Share on other sites

This will do:

internal void TestMath() {
	constant F32 allowedError = 0.00001f;
	constant F32 cs45 = 0.707106769f; // NOTE(finaL): 45 degrees unit axis scalar

	F32 sResult;

	sResult = DegreeToRadians32(45);
	Assert(Equals(sResult, PI32 * 0.25f, allowedError));
	sResult = RadiansToDegree32(PI32 * 0.25f);
	Assert(Equals(sResult, 45.0f, allowedError));

	sResult = ToAngle(V2(0, 1));
	Assert(Equals(sResult, PI32 * 0.5f, allowedError));
	sResult = ToAngle(V2(0, -1));
	Assert(Equals(sResult, PI32 * -0.5f, allowedError));
	sResult = ToAngle(V2(1, 0));
	Assert(Equals(sResult, 0, allowedError));
	sResult = ToAngle(V2(-1, 0));
	Assert(Equals(sResult, PI32, allowedError));
	sResult = ToAngle(V2(cs45, cs45));
	Assert(Equals(sResult, PI32 * 0.25f, allowedError));

	Mat2f rotIdentity = MakeIdentity2();
	Assert(
		Equals(rotIdentity.col1.x, 1, allowedError)
		&& Equals(rotIdentity.col1.y, 0, allowedError)
		&& Equals(rotIdentity.col2.x, 0, allowedError)
		&& Equals(rotIdentity.col2.y, 1, allowedError));

	Mat2f rot90 = MakeRotation2(PI32 * 0.5f);
	Assert(
		Equals(rot90.col1.x, 0, allowedError)
		&& Equals(rot90.col1.y, 1, allowedError)
		&& Equals(rot90.col2.x, -1, allowedError)
		&& Equals(rot90.col2.y, 0, allowedError));

	Mat2f rot45 = MakeRotation2(PI32 * 0.25f);
	Assert(
		Equals(rot45.col1.x, 0.707106769f, allowedError)
		&& Equals(rot45.col1.y, 0.707106769f, allowedError)
		&& Equals(rot45.col2.x, -0.707106769f, allowedError)
		&& Equals(rot45.col2.y, 0.707106769f, allowedError));

	Mat2f rot90Inverse = Transpose2(rot90);
	Assert(
		Equals(rot90Inverse.col1.x, 0, allowedError)
		&& Equals(rot90Inverse.col1.y, -1, allowedError)
		&& Equals(rot90Inverse.col2.x, 1, allowedError)
		&& Equals(rot90Inverse.col2.y, 0, allowedError));

	Vec2f v, vResult;
	Mat2f mResult;

	v = V2(100, 25);

	// NOTE(final): Matrix * vector is equal to vector * matrix
	vResult = rotIdentity * v;
	Assert(Equals(vResult.x, 100.0f, allowedError) && Equals(vResult.y, 25.0f, allowedError));
	vResult = v * rotIdentity;
	Assert(Equals(vResult.x, 100.0f, allowedError) && Equals(vResult.y, 25.0f, allowedError));

	vResult = rot90 * v;
	Assert(Equals(vResult.x, -25.0f, allowedError) && Equals(vResult.y, 100.0f, allowedError));
	vResult = v * rot90;
	Assert(Equals(vResult.x, -25.0f, allowedError) && Equals(vResult.y, 100.0f, allowedError));

	vResult = rot90Inverse * v;
	Assert(Equals(vResult.x, 25.0f, allowedError) && Equals(vResult.y, -100.0f, allowedError));
	vResult = v * rot90Inverse;
	Assert(Equals(vResult.x, 25.0f, allowedError) && Equals(vResult.y, -100.0f, allowedError));

	mResult = rot90Inverse * rot90;
	Assert(
		Equals(mResult.col1.x, 1, allowedError)
		&& Equals(mResult.col1.y, 0, allowedError)
		&& Equals(mResult.col2.x, 0, allowedError)
		&& Equals(mResult.col2.y, 1, allowedError));

	mResult = rot90 * rot90Inverse;
	Assert(
		Equals(mResult.col1.x, 1, allowedError)
		&& Equals(mResult.col1.y, 0, allowedError)
		&& Equals(mResult.col2.x, 0, allowedError)
		&& Equals(mResult.col2.y, 1, allowedError));

	mResult = rot45 * rot90;
	Assert(
		Equals(mResult.col1.x, -cs45, allowedError)
		&& Equals(mResult.col1.y, cs45, allowedError)
		&& Equals(mResult.col2.x, -cs45, allowedError)
		&& Equals(mResult.col2.y, -cs45, allowedError));
}

Share this post


Link to post
Share on other sites

This topic is 416 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