/Utilities/Datatypes/Matrix.cs

using System;

using System.IO;

using System.Runtime.InteropServices;

using Delta.Utilities.Helpers;

using Delta.Utilities.Profiling;

using NUnit.Framework;





// Disable warnings for some tests where we don't use created values

#pragma warning disable 219

// Disable the obsolete warnings for Matrix.Multiply here, it is optimized

// already here and when it is not in some unit test, its ok too.

#pragma warning disable 618



namespace Delta.Utilities.Datatypes

{

	/// <summary>

	/// Matrix struct, used mostly for 3D calculations. If possible this shares

	/// all the functionality of the XNA Matrix struct (if we are compiling

	/// only). In case we do not compile with XNA or SlimDX enabled (like on

	/// this class still provides all the basic functionality as a fallback. It

	/// might not be as optimized as some specialized xna code paths for the

	/// Xbox 360, but those only work on the XNA define and platform anyway.

	/// </summary>

	/// <remarks>

	/// Represents a row based matrix. Good wiki page about the coordinate

	/// system (right handed):

	/// http://en.wikipedia.org/wiki/Cartesian_coordinate_system

	/// <para />

	/// Note: This Matrix class was heavily profiled and optimized with help of

	/// the MatrixPerformance helper class included here. Many code paths

	/// are platform and framework dependant, this is why this class is so huge,

	/// but since this is very low level and highly performance critical, it is

	/// well worth the complexity and effort (the class is still easy to use).

	/// On the Xbox we will always try to use XNA implementation for performance

	/// critical methods, all other platforms are mostly implemented by us

	/// after testing performance. You can also checkout the SlimDX math classes,

	/// which mostly are just C#, only certain heavy methods like multiplying

	/// an array of matrices uses PInvoke to D3DXMatrixMultiply for example.

	/// <para />

	/// There are still performance improvements possible for this class

	/// like using ref for most performance critical methods like Invert,

	/// Multiply and most other operators, but that makes this class much harder

	/// to use and we can still do code-optimizations directly with the exposed

	/// data when needed (e.g. Bone Matrix code should be heavily optimized if

	/// it has to run on CPU). Currently this is fixed by marking slow methods

	/// as obsolete, which will give you compiler warnings, use the ref

	/// overload methods instead.

	/// <para />

	/// Also note that many methods are not faster by using the native XNA,

	/// OpenTK or SharpDX or SlimDX code paths, especially when just assigning

	/// values (where our code is almost always faster). But remember when any

	/// framework is using special tricks like low level assembly code or native

	/// code to calculate complex math operations, it can be faster and will be

	/// replaced by the build system automatically. If you notice any method

	/// that should be faster, contact us immediately so we can investigate!

	/// </remarks>

	[Serializable]

	[StructLayout(LayoutKind.Explicit)]

	public struct Matrix : ISaveLoadBinary, IEquatable<Matrix>

	{

		#region Constants

		/// <summary>

		/// Represents the number of values of that Matrix struct (here 4x4 = 16).

		/// </summary>

		public const int ValueCount = 16;



		/// <summary>

		/// Returns a zeroed matrix.

		/// </summary>

		public static readonly Matrix Zero = new Matrix(

			0, 0, 0, 0,

			0, 0, 0, 0,

			0, 0, 0, 0,

			0, 0, 0, 0);

		#endregion



		#region Invert (Static)

		/// <summary>

		/// Invert matrix, this is only slightly faster on the xna platform

		/// specific implementation, for performance critical code use the ref

		/// version!

		/// </summary>

		/// <param name="a">Matrix to invert</param>

		/// <returns>Inverted Matrix</returns>

		public static Matrix Invert(Matrix a)

		{

			//Check performance difference again:



			//Note: Performance results from 10 million calls:

			//Delta: 1962ms

			//XNA: 1076ms

			//OpenTK: 8611ms

			//SlimDX: 981ms



			// Ugly code that is actually lots faster than our clean implementation

			// below (see link on how it works basically). This is basically getting

			// the Determinant and then also calculating the inverse matrix in one

			// big ugly swoop.

			float m11 = a.M11;

			float m12 = a.M12;

			float m13 = a.M13;

			float m14 = a.M14;

			float m21 = a.M21;

			float m22 = a.M22;

			float m23 = a.M23;

			float m24 = a.M24;

			float m31 = a.M31;

			float m32 = a.M32;

			float m33 = a.M33;

			float m34 = a.M34;

			float m41 = a.M41;

			float m42 = a.M42;

			float m43 = a.M43;

			float m44 = a.M44;

			float m33x44m34x43 = (m33 * m44) - (m34 * m43);

			float m32x44m34x42 = (m32 * m44) - (m34 * m42);

			float m32x43m33x42 = (m32 * m43) - (m33 * m42);

			float m31x44m34x41 = (m31 * m44) - (m34 * m41);

			float m31x43m33x41 = (m31 * m43) - (m33 * m41);

			float m31x42m32x41 = (m31 * m42) - (m32 * m41);

			float m22m23m24 = ((m22 * m33x44m34x43) - (m23 * m32x44m34x42)) +

			                  (m24 * m32x43m33x42);

			float m21m23m24 = -(((m21 * m33x44m34x43) - (m23 * m31x44m34x41)) +

			                    (m24 * m31x43m33x41));

			float m21m22m24 = ((m21 * m32x44m34x42) - (m22 * m31x44m34x41)) +

			                  (m24 * m31x42m32x41);

			float m21m22m23 = -(((m21 * m32x43m33x42) - (m22 * m31x43m33x41)) +

			                    (m23 * m31x42m32x41));



			// Compute the the determinant

			float det =

				(m11 * m22m23m24) + (m12 * m21m23m24) + (m13 * m21m22m24) +

				(m14 * m21m22m23);

			if (det == 0.0f)

			{

				det = 0.000001f;

			}



			// Inverse(A) = 1/det(A) * B

			float inverseDet = 1f / det;

			float m23x44m24x43 = (m23 * m44) - (m24 * m43);

			float m22x44m24x42 = (m22 * m44) - (m24 * m42);

			float m22x43m23x42 = (m22 * m43) - (m23 * m42);

			float m21x44m24x41 = (m21 * m44) - (m24 * m41);

			float m21x43m23x41 = (m21 * m43) - (m23 * m41);

			float m21x42m22x41 = (m21 * m42) - (m22 * m41);

			float m23x34m24x33 = (m23 * m34) - (m24 * m33);

			float m22x34m24x32 = (m22 * m34) - (m24 * m32);

			float m22x33m23x32 = (m22 * m33) - (m23 * m32);

			float m21x34m23x31 = (m21 * m34) - (m24 * m31);

			float m21x33m23x31 = (m21 * m33) - (m23 * m31);

			float m21x32m22x31 = (m21 * m32) - (m22 * m31);

			return new Matrix(

				m22m23m24 * inverseDet,

				-(((m12 * m33x44m34x43) - (m13 * m32x44m34x42)) +

				  (m14 * m32x43m33x42)) * inverseDet,

				(((m12 * m23x44m24x43) - (m13 * m22x44m24x42)) +

				 (m14 * m22x43m23x42)) * inverseDet,

				-(((m12 * m23x34m24x33) - (m13 * m22x34m24x32)) +

				  (m14 * m22x33m23x32)) * inverseDet,

				m21m23m24 * inverseDet,

				(((m11 * m33x44m34x43) - (m13 * m31x44m34x41)) +

				 (m14 * m31x43m33x41)) * inverseDet,

				-(((m11 * m23x44m24x43) - (m13 * m21x44m24x41)) +

				  (m14 * m21x43m23x41)) * inverseDet,

				(((m11 * m23x34m24x33) - (m13 * m21x34m23x31)) +

				 (m14 * m21x33m23x31)) * inverseDet,

				m21m22m24 * inverseDet,

				-(((m11 * m32x44m34x42) - (m12 * m31x44m34x41)) +

				  (m14 * m31x42m32x41)) * inverseDet,

				(((m11 * m22x44m24x42) - (m12 * m21x44m24x41)) +

				 (m14 * m21x42m22x41)) * inverseDet,

				-(((m11 * m22x34m24x32) - (m12 * m21x34m23x31)) +

				  (m14 * m21x32m22x31)) * inverseDet,

				m21m22m23 * inverseDet,

				(((m11 * m32x43m33x42) - (m12 * m31x43m33x41)) +

				 (m13 * m31x42m32x41)) * inverseDet,

				-(((m11 * m22x43m23x42) - (m12 * m21x43m23x41)) +

				  (m13 * m21x42m22x41)) * inverseDet,

				(((m11 * m22x33m23x32) - (m12 * m21x33m23x31)) +

				 (m13 * m21x32m22x31)) * inverseDet);

		}



		/// <summary>

		/// Invert matrix ref version, this is the fastest one, especially if

		/// you use matrix and result for the same value, which will copy values

		/// over quickly.

		/// </summary>

		/// <param name="matrix">Matrix</param>

		/// <param name="result">Result</param>

		public static void Invert(ref Matrix matrix, ref Matrix result)

		{

			//Check performance difference again:



			// Ugly code that is actually lots faster than our clean implementation

			// below (see link on how it works basically). This is basically getting

			// the Determinant and then also calculating the inverse matrix in one

			// big ugly swoop.

			float m11 = matrix.M11;

			float m12 = matrix.M12;

			float m13 = matrix.M13;

			float m14 = matrix.M14;

			float m21 = matrix.M21;

			float m22 = matrix.M22;

			float m23 = matrix.M23;

			float m24 = matrix.M24;

			float m31 = matrix.M31;

			float m32 = matrix.M32;

			float m33 = matrix.M33;

			float m34 = matrix.M34;

			float m41 = matrix.M41;

			float m42 = matrix.M42;

			float m43 = matrix.M43;

			float m44 = matrix.M44;

			float m33x44m34x43 = (m33 * m44) - (m34 * m43);

			float m32x44m34x42 = (m32 * m44) - (m34 * m42);

			float m32x43m33x42 = (m32 * m43) - (m33 * m42);

			float m31x44m34x41 = (m31 * m44) - (m34 * m41);

			float m31x43m33x41 = (m31 * m43) - (m33 * m41);

			float m31x42m32x41 = (m31 * m42) - (m32 * m41);

			float m22m23m24 = ((m22 * m33x44m34x43) - (m23 * m32x44m34x42)) +

			                  (m24 * m32x43m33x42);

			float m21m23m24 = -(((m21 * m33x44m34x43) - (m23 * m31x44m34x41)) +

			                    (m24 * m31x43m33x41));

			float m21m22m24 = ((m21 * m32x44m34x42) - (m22 * m31x44m34x41)) +

			                  (m24 * m31x42m32x41);

			float m21m22m23 = -(((m21 * m32x43m33x42) - (m22 * m31x43m33x41)) +

			                    (m23 * m31x42m32x41));



			// Computing the determinant

			float det = (m11 * m22m23m24) + (m12 * m21m23m24) + (m13 * m21m22m24) +

			            (m14 * m21m22m23);

			if (det == 0.0f)

			{

				det = 0.000001f;

			}



			// Inverse(A) = 1/det(A) * B

			float inverseDet = 1f / det;

			float m23x44m24x43 = (m23 * m44) - (m24 * m43);

			float m22x44m24x42 = (m22 * m44) - (m24 * m42);

			float m22x43m23x42 = (m22 * m43) - (m23 * m42);

			float m21x44m24x41 = (m21 * m44) - (m24 * m41);

			float m21x43m23x41 = (m21 * m43) - (m23 * m41);

			float m21x42m22x41 = (m21 * m42) - (m22 * m41);

			float m23x34m24x33 = (m23 * m34) - (m24 * m33);

			float m22x34m24x32 = (m22 * m34) - (m24 * m32);

			float m22x33m23x32 = (m22 * m33) - (m23 * m32);

			float m21x34m23x31 = (m21 * m34) - (m24 * m31);

			float m21x33m23x31 = (m21 * m33) - (m23 * m31);

			float m21x32m22x31 = (m21 * m32) - (m22 * m31);

			// Now just set everything to the output result matrix, which even can

			// be the same as the input matrix (for best performance).

			result.M11 = m22m23m24 * inverseDet;

			result.M12 = -(((m12 * m33x44m34x43) - (m13 * m32x44m34x42)) +

			               (m14 * m32x43m33x42)) * inverseDet;

			result.M13 = (((m12 * m23x44m24x43) - (m13 * m22x44m24x42)) +

			              (m14 * m22x43m23x42)) * inverseDet;

			result.M14 = -(((m12 * m23x34m24x33) - (m13 * m22x34m24x32)) +

			               (m14 * m22x33m23x32)) * inverseDet;

			result.M21 = m21m23m24 * inverseDet;

			result.M22 = (((m11 * m33x44m34x43) - (m13 * m31x44m34x41)) +

			              (m14 * m31x43m33x41)) * inverseDet;

			result.M23 = -(((m11 * m23x44m24x43) - (m13 * m21x44m24x41)) +

			               (m14 * m21x43m23x41)) * inverseDet;

			result.M24 = (((m11 * m23x34m24x33) - (m13 * m21x34m23x31)) +

			              (m14 * m21x33m23x31)) * inverseDet;

			result.M31 = m21m22m24 * inverseDet;

			result.M32 = -(((m11 * m32x44m34x42) - (m12 * m31x44m34x41)) +

			               (m14 * m31x42m32x41)) * inverseDet;

			result.M33 = (((m11 * m22x44m24x42) - (m12 * m21x44m24x41)) +

			              (m14 * m21x42m22x41)) * inverseDet;

			result.M34 = -(((m11 * m22x34m24x32) - (m12 * m21x34m23x31)) +

			               (m14 * m21x32m22x31)) * inverseDet;

			result.M41 = m21m22m23 * inverseDet;

			result.M42 = (((m11 * m32x43m33x42) - (m12 * m31x43m33x41)) +

			              (m13 * m31x42m32x41)) * inverseDet;

			result.M43 = -(((m11 * m22x43m23x42) - (m12 * m21x43m23x41)) +

			               (m13 * m21x42m22x41)) * inverseDet;

			result.M44 = (((m11 * m22x33m23x32) - (m12 * m21x33m23x31)) +

			              (m13 * m21x32m22x31)) * inverseDet;

		}

		#endregion



		#region Transpose (Static)

		/// <summary>

		/// Transpose, which basically just means we switch from a row to a column

		/// based matrix.

		/// </summary>

		/// <param name="matrix">Matrix</param>

		/// <returns>Transposed matrix (row/columns switched)</returns>

		public static Matrix Transpose(Matrix matrix)

		{

			//Check performance difference again:



			// Note: In all cases our implementation is slighly faster/or a lot faster

			//Delta: 569ms (10 mio calls)

			//XNA: 760ms (10 mio calls)

			//OpenTK: 904ms (10 mio calls)

			//SlimDX: 658ms (10 mio calls)



			return new Matrix(

				// First row

				matrix.M11, matrix.M21, matrix.M31, matrix.M41,

				// Second row

				matrix.M12, matrix.M22, matrix.M32, matrix.M42,

				// Third row

				matrix.M13, matrix.M23, matrix.M33, matrix.M43,

				// Forth row

				matrix.M14, matrix.M24, matrix.M34, matrix.M44);

		}



		/// <summary>

		/// Transpose, which basically just means we switch from a row to a column

		/// based matrix. Faster version using ref input and output values.

		/// </summary>

		/// <param name="matrix">Matrix</param>

		/// <param name="result">result</param>

		public static void Transpose(ref Matrix matrix, ref Matrix result)

		{

			result.M11 = matrix.M11;

			result.M12 = matrix.M21;

			result.M13 = matrix.M31;

			result.M14 = matrix.M41;

			result.M21 = matrix.M12;

			result.M22 = matrix.M22;

			result.M23 = matrix.M32;

			result.M24 = matrix.M42;

			result.M31 = matrix.M13;

			result.M32 = matrix.M23;

			result.M33 = matrix.M33;

			result.M34 = matrix.M43;

			result.M41 = matrix.M14;

			result.M42 = matrix.M24;

			result.M43 = matrix.M34;

			result.M44 = matrix.M44;

		}

		#endregion



		#region CreateColumnBased (Static)

		/// <summary>

		/// Create column based matrix, also used for creating an exported collada

		/// matrix, because the normal XNA/OpenTK/whatever matrix is row based.

		/// </summary>

		/// <param name="floatValues">Float values</param>

		/// <returns>Matrix</returns>

		public static Matrix CreateColumnBased(float[] floatValues)

		{

			// If we get no valid or an incorrect number of float values, then we

			// return only the identity matrix

			if (floatValues == null ||

			    floatValues.Length != 16)

			{

				return Identity;

			}



			return new Matrix(

				floatValues[0], floatValues[4], floatValues[8], floatValues[12],

				floatValues[1], floatValues[5], floatValues[9], floatValues[13],

				floatValues[2], floatValues[6], floatValues[10], floatValues[14],

				floatValues[3], floatValues[7], floatValues[11], floatValues[15]);

		}

		#endregion



		#region CreateScale (Static)

		/// <summary>

		/// Create a scale matrix with the specified value for X, Y and Z.

		/// </summary>

		/// <param name="scale">Scale</param>

		/// <returns>New Matrix with the given scale</returns>

		public static Matrix CreateScale(float scale)

		{

			return new Matrix(

				scale, 0f, 0f, 0f,

				0f, scale, 0f, 0f,

				0f, 0f, scale, 0f,

				0f, 0f, 0f, 1.0f);

		}



		/// <summary>

		/// Create scale

		/// </summary>

		/// <param name="scale">Scale</param>

		/// <returns>New Matrix with the given scale</returns>

		public static Matrix CreateScale(Vector scale)

		{

			return CreateScale(scale.X, scale.Y, scale.Z);

		}



		/// <summary>

		/// Create a scale matrix with the specified X, Y and Z scaling values.

		/// </summary>

		/// <param name="scaleX">The amount of scaling in X dimension.</param>

		/// <param name="scaleY">The amount of scaling in Y dimension.</param>

		/// <param name="scaleZ">The amount of scaling in Z dimension.</param>

		/// <returns>New Matrix with the given scale</returns>

		public static Matrix CreateScale(float scaleX, float scaleY, float scaleZ)

		{

			return new Matrix(

				scaleX, 0f, 0f, 0f,

				0f, scaleY, 0f, 0f,

				0f, 0f, scaleZ, 0f,

				0f, 0f, 0f, 1f);

		}

		#endregion



		#region CreateRotationX (Static)

		/// <summary>

		/// Create rotation around x-axis (pitch)

		/// <para />

		/// It is important to know which axis is meant with Yaw, Pitch and Roll:

		/// http://www.spotimage.com/dimap/spec/dictionary/Spot_Scene/Illustrations/YAW.gif

		/// </summary>

		public static Matrix CreateRotationX(float degrees)

		{

			//Check performance difference again:



			// Note: All CreateRotation implementations are slow, there are slightly

			// faster on XNA, but overall not much of a difference between all

			// platforms and frameworks (see below in MatrixPerformance class).



			float cosValue = MathHelper.Cos(degrees);

			float sinValue = MathHelper.Sin(degrees);

			return new Matrix(

				1f, 0f, 0f, 0f,

				0f, cosValue, sinValue, 0f,

				0f, -sinValue, cosValue, 0f,

				0f, 0f, 0f, 1f);

		}

		#endregion



		#region CreateRotationY (Static)

		/// <summary>

		/// Create rotation around y-axis (yaw)

		/// <para />

		/// It is important to know which axis is meant with Yaw, Pitch and Roll:

		/// http://www.spotimage.com/dimap/spec/dictionary/Spot_Scene/Illustrations/YAW.gif

		/// </summary>

		/// <param name="degrees">Degrees to rotate around the Y axis</param>

		/// <returns>New matrix with the given rotation</returns>

		public static Matrix CreateRotationY(float degrees)

		{

			//Check performance difference again:

			// Note: All CreateRotation implementations are slow, there are slightly

			// faster on XNA, but overall not much of a difference between all

			// platforms and frameworks (see below in MatrixPerformance class).



			float cosValue = MathHelper.Cos(degrees);

			float sinValue = MathHelper.Sin(degrees);

			return new Matrix(

				cosValue, 0f, -sinValue, 0f,

				0f, 1f, 0f, 0f,

				sinValue, 0f, cosValue, 0f,

				0f, 0f, 0f, 1f);

		}

		#endregion



		#region CreateRotationZ (Static)

		/// <summary>

		/// Create rotation around z-axis (roll)

		/// <para />

		/// It is important to know which axis is meant with Yaw, Pitch and Roll:

		/// http://www.spotimage.com/dimap/spec/dictionary/Spot_Scene/Illustrations/YAW.gif

		/// </summary>

		/// <param name="degrees">Degrees to rotate around the Z axis</param>

		/// <returns>New matrix with the given rotation</returns>

		public static Matrix CreateRotationZ(float degrees)

		{

			//Check performance difference again:



			// Note: All CreateRotation implementations are slow, there are slightly

			// faster on XNA, but overall not much of a difference between all

			// platforms and frameworks (see below in MatrixPerformance class).



			float cosValue = MathHelper.Cos(degrees);

			float sinValue = MathHelper.Sin(degrees);

			return new Matrix(

				cosValue, sinValue, 0f, 0f,

				-sinValue, cosValue, 0f, 0f,

				0f, 0f, 1f, 0f,

				0f, 0f, 0f, 1f);

		}

		#endregion



		#region CreateRotationZYX (Static)

		/// <summary>

		/// Rotate around the Z axis, then the Y axis, and finally the X axis

		/// (rotX * rotY * rotZ). Please note that this is not the same as

		/// FromYawPitchRoll, which uses Z (yaw), then X (pitch), then Y (roll).

		/// <para />

		/// Note: This is the same as calling the CreateRotationX, Y, Z methods,

		/// just combined all together without having to multiply matrices).

		/// </summary>

		/// <param name="x">Degrees to rotate around the X axis</param>

		/// <param name="y">Degrees to rotate around the Y axis</param>

		/// <param name="z">Degrees to rotate around the Z axis</param>

		/// <returns>New matrix with the given rotation</returns>

		public static Matrix CreateRotationZYX(float x, float y, float z)

		{

			float cx = 1;

			float sx = 0;

			float cy = 1;

			float sy = 0;

			float cz = 1;

			float sz = 0;

			if (x != 0.0f)

			{

				cx = MathHelper.Cos(x);

				sx = MathHelper.Sin(x);

			}

			if (y != 0.0f)

			{

				cy = MathHelper.Cos(y);

				sy = MathHelper.Sin(y);

			}

			if (z != 0.0f)

			{

				cz = MathHelper.Cos(z);

				sz = MathHelper.Sin(z);

			}



			return new Matrix(

				cy * cz, cy * sz, -sy, 0.0f,

				(sx * sy * cz) + (cx * -sz), (sx * sy * sz) + (cx * cz), sx * cy, 0.0f,

				(cx * sy * cz) + (sx * sz), (cx * sy * sz) + (-sx * cz), cx * cy, 0.0f,

				0.0f, 0.0f, 0.0f, 1.0f);

		}

		#endregion



		#region CreateRotationZY (Static)

		/// <summary>

		/// Rotate around Z axis, then around Y axis.

		/// (equals rotY * rotZ)

		/// </summary>

		/// <param name="y">Degrees to rotate around the Y axis</param>

		/// <param name="z">Degrees to rotate around the Z axis</param>

		/// <returns>New matrix with the given rotation</returns>

		public static Matrix CreateRotationZY(float y, float z)

		{

			float cy = 1;

			float sy = 0;

			float cz = 1;

			float sz = 0;

			// We only calculate sin/cos if given input is not zero

			if (y != 0.0f)

			{

				cy = MathHelper.Cos(y);

				sy = MathHelper.Sin(y);

			}

			if (z != 0.0f)

			{

				cz = MathHelper.Cos(z);

				sz = MathHelper.Sin(z);

			}

			return new Matrix(

				cy * cz, cy * sz, -sy, 0.0f,

				-sz, cz, 0.0f, 0.0f,

				sy * cz, sy * sz, cy, 0.0f,

				0.0f, 0.0f, 0.0f, 1.0f);

		}

		#endregion



		#region FromYawPitchRoll (Static)

		/// <summary>

		/// From yaw pitch roll (yaw = Y, pitch = X, roll = Z)

		/// http://www.euclideanspace.com/maths/geometry/rotations/conversions/eulerToMatrix/index.htm

		/// It is important to know which axis is meant with Yaw, Pitch and Roll:

		/// http://www.spotimage.com/dimap/spec/dictionary/Spot_Scene/Illustrations/YAW.gif

		/// <para />

		/// Note: This method calculates Y * X * Z where

		/// Y = value.X = yaw

		/// X = value.Y = pitch

		/// Z = value.Z = roll

		/// </summary>

		public static Matrix FromYawPitchRoll(Vector value)

		{

			return FromYawPitchRoll(value.X, value.Y, value.Z);

		}



		/// <summary>

		/// From yaw pitch roll (yaw = Y, pitch = X, roll = Z)

		/// It is important to know which axis is meant with Yaw, Pitch and Roll:

		/// http://www.spotimage.com/dimap/spec/dictionary/Spot_Scene/Illustrations/YAW.gif

		/// <para />

		/// Note: This method calculates Y * X * Z where

		/// Y = yaw

		/// X = pitch

		/// Z = roll

		/// </summary>

		public static Matrix FromYawPitchRoll(float yaw, float pitch, float roll)

		{

			//Check performance difference again:



			float cy = MathHelper.Cos(yaw);

			float sy = MathHelper.Sin(yaw);

			float cp = MathHelper.Cos(pitch);

			float sp = MathHelper.Sin(pitch);

			float cr = MathHelper.Cos(roll);

			float sr = MathHelper.Sin(roll);



			return new Matrix(

				cy * cr - sr * sy * sp, cy * sr + sy * sp * cr, -sy * cp, 0.0f,

				-sr * cp, cp * cr, sp, 0.0f,

				sy * cr + sp * cy * sr, sy * sr - sp * cy * cr, cy * cp, 0.0f,

				0.0f, 0.0f, 0.0f, 1.0f);

		}

		#endregion



		#region CreateTranslation (Static)

		/// <summary>

		/// Create translation matrix

		/// </summary>

		/// <param name="position">Position</param>

		/// <returns>Identity matrix with the translation (M41, M42, M43)</returns>

		public static Matrix CreateTranslation(Vector position)

		{

			return new Matrix(

				1f, 0f, 0f, 0f,

				0f, 1f, 0f, 0f,

				0f, 0f, 1f, 0f,

				position.X, position.Y, position.Z, 1.0f);

		}



		/// <summary>

		/// Create translation matrix

		/// </summary>

		/// <param name="x">x</param>

		/// <param name="y">y</param>

		/// <param name="z">z</param>

		/// <returns>Identity matrix with the translation (M41, M42, M43)</returns>

		public static Matrix CreateTranslation(float x, float y, float z)

		{

			return new Matrix(

				1f, 0f, 0f, 0f,

				0f, 1f, 0f, 0f,

				0f, 0f, 1f, 0f,

				x, y, z, 1.0f);

		}

		#endregion



		#region CreateScaleAndTranslation (Static)

		/// <summary>

		/// Create scale and translation, basically just does CreateScale and

		/// CreateTranslation, but much more efficient.

		/// </summary>

		/// <param name="scale">Scale for the new matrix</param>

		/// <param name="position">Position for the matrix translation.</param>

		/// <returns>New matrix with the given scale and translation.</returns>

		public static Matrix CreateScaleAndTranslation(float scale,

			Vector position)

		{

			return new Matrix(

				scale, 0f, 0f, 0f,

				0f, scale, 0f, 0f,

				0f, 0f, scale, 0f,

				position.X, position.Y, position.Z, 1f);

		}

		#endregion



		#region CreateLookAt (Static)

		/// <summary>

		/// Creates a right-handed look at matrix for a camera.

		/// </summary>

		/// <param name="cameraPosition">camera position</param>

		/// <param name="cameraTarget">camera target</param>

		/// <param name="cameraUpVector">camera up vector</param>

		/// <returns>New look at matrix</returns>

		public static Matrix CreateLookAt(Vector cameraPosition,

			Vector cameraTarget, Vector cameraUpVector)

		{

			//Check performance difference again:



			// Slightly faster in XNA (5%), but overall all these implementations

			// are roughly the same.

			// Note: All Create view implementations are slow, there are slighly

			// faster on XNA, but overall not much of a difference between all

			// platforms and frameworks (see below in MatrixPerformance class).





			// Simplified and optimized formula to create look matrix

			Vector dir = Vector.Normalize(cameraPosition - cameraTarget);

			Vector up = Vector.Normalize(Vector.Cross(cameraUpVector, dir));

			// Switch up around if cameraUpVector is the same as dir!

			if (up.LengthSquared == 0)

			{

				up = Vector.UnitY;

			}

			Vector right = Vector.Cross(dir, up);

			return new Matrix(

				up.X, right.X, dir.X, 0f,

				up.Y, right.Y, dir.Y, 0f,

				up.Z, right.Z, dir.Z, 0f,

				-Vector.Dot(up, cameraPosition),

				-Vector.Dot(right, cameraPosition),

				-Vector.Dot(dir, cameraPosition), 1f);

		}

		#endregion



		#region CreateFromAxisAngle (Static)

		/// <summary>

		/// Create from axis angle (in degrees)

		/// </summary>

		/// <param name="angle">Angle in degrees</param>

		/// <param name="axis">Axis to rotate around</param>

		/// <returns>Rotated matrix around axis for rotations</returns>

		public static Matrix CreateFromAxisAngle(Vector axis, float angle)

		{

			// Note: Variables in this method need to be beautifulized!

			//Check performance difference again:



			Matrix matrix = new Matrix();



			float x = axis.X;

			float y = axis.Y;

			float z = axis.Z;

			float num2 = MathHelper.Sin(angle);

			float num = MathHelper.Cos(angle);

			float num11 = x * x;

			float num10 = y * y;

			float num9 = z * z;

			float num8 = x * y;

			float num7 = x * z;

			float num6 = y * z;



			matrix.M11 = num11 + (num * (1f - num11));

			matrix.M12 = (num8 - (num * num8)) + (num2 * z);

			matrix.M13 = (num7 - (num * num7)) - (num2 * y);

			matrix.M14 = 0f;

			matrix.M21 = (num8 - (num * num8)) - (num2 * z);

			matrix.M22 = num10 + (num * (1f - num10));

			matrix.M23 = (num6 - (num * num6)) + (num2 * x);

			matrix.M24 = 0f;

			matrix.M31 = (num7 - (num * num7)) + (num2 * y);

			matrix.M32 = (num6 - (num * num6)) - (num2 * x);

			matrix.M33 = num9 + (num * (1f - num9));

			matrix.M34 = 0f;

			matrix.M41 = 0f;

			matrix.M42 = 0f;

			matrix.M43 = 0f;

			matrix.M44 = 1f;



			return matrix;

		}



		/// <summary>

		/// Create from axis angle (in degrees)

		/// </summary>

		/// <param name="angle">Angle in degrees</param>

		/// <param name="axis">Axis to rotate around</param>

		/// <param name="result">Rotated matrix around axis for rotations</param>

		public static void CreateFromAxisAngle(ref Vector axis, float angle,

			ref Matrix result)

		{

			float x = axis.X;

			float y = axis.Y;

			float z = axis.Z;

			float num2 = MathHelper.Sin(angle);

			float num = MathHelper.Cos(angle);

			float num11 = x * x;

			float num10 = y * y;

			float num9 = z * z;

			float num8 = x * y;

			float num7 = x * z;

			float num6 = y * z;



			result.M11 = num11 + (num * (1f - num11));

			result.M12 = (num8 - (num * num8)) + (num2 * z);

			result.M13 = (num7 - (num * num7)) - (num2 * y);

			result.M14 = 0f;

			result.M21 = (num8 - (num * num8)) - (num2 * z);

			result.M22 = num10 + (num * (1f - num10));

			result.M23 = (num6 - (num * num6)) + (num2 * x);

			result.M24 = 0f;

			result.M31 = (num7 - (num * num7)) + (num2 * y);

			result.M32 = (num6 - (num * num6)) - (num2 * x);

			result.M33 = num9 + (num * (1f - num9));

			result.M34 = 0f;

			result.M41 = 0f;

			result.M42 = 0f;

			result.M43 = 0f;

			result.M44 = 1f;

		}

		#endregion



		#region CreateOrthographic (Static)

		/// <summary>

		/// Create orthographic

		/// </summary>

		/// <param name="bottom">bottom</param>

		/// <param name="farPlane">farPlane</param>

		/// <param name="left">left</param>

		/// <param name="nearPlane">near plane</param>

		/// <param name="right">right</param>

		/// <param name="top">top</param>

		/// <returns>New Orthographic matrix</returns>

		public static Matrix CreateOrthographic(float left, float right,

			float bottom, float top, float nearPlane, float farPlane)

		{

			//Check performance difference again:



			Matrix matrix = new Matrix();



			matrix.M11 = 2f / (right - left);

			matrix.M12 = matrix.M13 = matrix.M14 = 0f;

			matrix.M22 = 2f / (top - bottom);

			matrix.M21 = matrix.M23 = matrix.M24 = 0f;

			matrix.M33 = 1f / (nearPlane - farPlane);

			matrix.M31 = matrix.M32 = matrix.M34 = 0f;

			matrix.M41 = (left + right) / (left - right);

			matrix.M42 = (top + bottom) / (bottom - top);

			matrix.M43 = nearPlane / (nearPlane - farPlane);

			matrix.M44 = 1f;



			return matrix;

		}

		#endregion



		#region CreatePerspective (Static)

		/// <summary>

		/// Creates a right-handed perspective field of view matrix for a camera.

		/// </summary>

		/// <param name="fieldOfView">Field of view (Warning: In radians)</param>

		/// <param name="aspectRatio">Aspect ratio</param>

		/// <param name="nearPlaneDistance">Near plane distance</param>

		/// <param name="farPlaneDistance">Far plane distance</param>

		/// <returns>Matrix</returns>

		public static Matrix CreatePerspective(float fieldOfView,

			float aspectRatio, float nearPlaneDistance, float farPlaneDistance)

		{

			//Check performance difference again:

			// Slightly faster in Delta (10%), but overall all this implementations

			// are ruffly the same. Same as the following, just faster: 

			//XnaMatrix.CreatePerspectiveFieldOfView(...)

			//TKMatrix.CreatePerspectiveFieldOfView(...)

			//SlimDX.Matrix.PerspectiveFovRH(...)



			float m22 = 1f / ((float)Math.Tan(fieldOfView * 0.5f));

			return new Matrix(

				m22 / aspectRatio, 0f, 0f, 0f,

				0f, m22, 0f, 0f,

				0f, 0f, farPlaneDistance / (nearPlaneDistance - farPlaneDistance), -1f,

				0f, 0f, (nearPlaneDistance * farPlaneDistance) /

				        (nearPlaneDistance - farPlaneDistance), 0f);

		}

		#endregion



		#region Multiply (Static)

		/// <summary>

		/// Ugly Matrix Multiply method (normally you would just use

		/// matrix1*matrix2 to multiply two matrices), but this one is a little

		/// faster because we don't need to copy over the matrix values and we can

		/// even provide a faster fallback on the xbox platform (using xna).

		/// Note: matrix1 or matrix2 can be the same as result, which results in

		/// even better performance for performance critical matrix multiply code.

		/// </summary>

		/// <param name="matrix1">matrix1</param>

		/// <param name="matrix2">matrix2</param>

		/// <param name="result">result</param>

		public static void Multiply(ref Matrix matrix1, ref Matrix matrix2,

			ref Matrix result)

		{

			//Check performance difference again:

			// Note: Values are stored in local variables to allow matrix1 or matrix2

			// and result be the same matrices (we won't overwrite M11-M44 while

			// calculating the result), which happens quite often in optimized code!



			float m11 = matrix1.M11 * matrix2.M11 + matrix1.M12 * matrix2.M21 +

			            matrix1.M13 * matrix2.M31 + matrix1.M14 * matrix2.M41;

			float m12 = matrix1.M11 * matrix2.M12 + matrix1.M12 * matrix2.M22 +

			            matrix1.M13 * matrix2.M32 + matrix1.M14 * matrix2.M42;

			float m13 = matrix1.M11 * matrix2.M13 + matrix1.M12 * matrix2.M23 +

			            matrix1.M13 * matrix2.M33 + matrix1.M14 * matrix2.M43;

			float m14 = matrix1.M11 * matrix2.M14 + matrix1.M12 * matrix2.M24 +

			            matrix1.M13 * matrix2.M34 + matrix1.M14 * matrix2.M44;

			float m21 = matrix1.M21 * matrix2.M11 + matrix1.M22 * matrix2.M21 +

			            matrix1.M23 * matrix2.M31 + matrix1.M24 * matrix2.M41;

			float m22 = matrix1.M21 * matrix2.M12 + matrix1.M22 * matrix2.M22 +

			            matrix1.M23 * matrix2.M32 + matrix1.M24 * matrix2.M42;

			float m23 = matrix1.M21 * matrix2.M13 + matrix1.M22 * matrix2.M23 +

			            matrix1.M23 * matrix2.M33 + matrix1.M24 * matrix2.M43;

			float m24 = matrix1.M21 * matrix2.M14 + matrix1.M22 * matrix2.M24 +

			            matrix1.M23 * matrix2.M34 + matrix1.M24 * matrix2.M44;

			float m31 = matrix1.M31 * matrix2.M11 + matrix1.M32 * matrix2.M21 +

			            matrix1.M33 * matrix2.M31 + matrix1.M34 * matrix2.M41;

			float m32 = matrix1.M31 * matrix2.M12 + matrix1.M32 * matrix2.M22 +

			            matrix1.M33 * matrix2.M32 + matrix1.M34 * matrix2.M42;

			float m33 = matrix1.M31 * matrix2.M13 + matrix1.M32 * matrix2.M23 +

			            matrix1.M33 * matrix2.M33 + matrix1.M34 * matrix2.M43;

			float m34 = matrix1.M31 * matrix2.M14 + matrix1.M32 * matrix2.M24 +

			            matrix1.M33 * matrix2.M34 + matrix1.M34 * matrix2.M44;

			float m41 = matrix1.M41 * matrix2.M11 + matrix1.M42 * matrix2.M21 +

			            matrix1.M43 * matrix2.M31 + matrix1.M44 * matrix2.M41;

			float m42 = matrix1.M41 * matrix2.M12 + matrix1.M42 * matrix2.M22 +

			            matrix1.M43 * matrix2.M32 + matrix1.M44 * matrix2.M42;

			float m43 = matrix1.M41 * matrix2.M13 + matrix1.M42 * matrix2.M23 +

			            matrix1.M43 * matrix2.M33 + matrix1.M44 * matrix2.M43;

			float m44 = matrix1.M41 * matrix2.M14 + matrix1.M42 * matrix2.M24 +

			            matrix1.M43 * matrix2.M34 + matrix1.M44 * matrix2.M44;

			result.M11 = m11;

			result.M12 = m12;

			result.M13 = m13;

			result.M14 = m14;

			result.M21 = m21;

			result.M22 = m22;

			result.M23 = m23;

			result.M24 = m24;

			result.M31 = m31;

			result.M32 = m32;

			result.M33 = m33;

			result.M34 = m34;

			result.M41 = m41;

			result.M42 = m42;

			result.M43 = m43;

			result.M44 = m44;

		}



		/// <summary>

		/// Ugly Matrix*Vector multiply method (basically transforming a vector

		/// with a matrix). Normally you would just use matrix*position to

		/// multiply, but this one is a little faster because we don't need to

		/// copy over the matrix values and we can even provide a faster fallback

		/// on the xbox platform (using xna).

		/// Note: position and result can be the same as result, which results in

		/// even better performance for performance critical matrix multiply code.

		/// </summary>

		/// <param name="matrix">matrix</param>

		/// <param name="position">position</param>

		/// <param name="result">result</param>

		public static void Multiply(ref Matrix matrix, ref Vector position,

			ref Vector result)

		{

			//Check performance difference again:



			float x = position.X * matrix.M11 + position.Y * matrix.M21 +

			          position.Z * matrix.M31 + matrix.M41;

			float y = position.X * matrix.M12 + position.Y * matrix.M22 +

			          position.Z * matrix.M32 + matrix.M42;

			float z = position.X * matrix.M13 + position.Y * matrix.M23 +

			          position.Z * matrix.M33 + matrix.M43;

			result.X = x;

			result.Y = y;

			result.Z = z;

		}



		/// <summary>

		/// Multiply matrix with a scale factor. Note: This is the slow version.

		/// This method needs to copy 128 bytes around (2 matrices, each 64 bytes).

		/// Try to use the ref version instead.

		/// </summary>

		/// <param name="matrix">Matrix</param>

		/// <param name="scaleFactor">Scale factor</param>

		/// <returns>

		/// New matrix with each value multiplied with scaleFactor

		/// </returns>

		public static Matrix Multiply(Matrix matrix, float scaleFactor)

		{

			return new Matrix(

				matrix.M11 * scaleFactor, matrix.M12 * scaleFactor,

				matrix.M13 * scaleFactor, matrix.M14 * scaleFactor,

				matrix.M21 * scaleFactor, matrix.M22 * scaleFactor,

				matrix.M23 * scaleFactor, matrix.M24 * scaleFactor,

				matrix.M31 * scaleFactor, matrix.M32 * scaleFactor,

				matrix.M33 * scaleFactor, matrix.M34 * scaleFactor,

				matrix.M41 * scaleFactor, matrix.M42 * scaleFactor,

				matrix.M43 * scaleFactor, matrix.M44 * scaleFactor);

		}

		#endregion



		#region Equal (Static)

		/// <summary>

		/// Equality check. Optimized for speed via ref parameters, which

		/// are much faster than copying 128 bytes (2 matrices each 64 bytes).

		/// </summary>

		public static bool Equal(ref Matrix matrix1, ref Matrix matrix2)

		{

			// Note: The checks are sorted by the diagonal, which is most important

			// and usually the most different for 2 matrices (faster this way!)

			// Note: Calling XNA (and OpenTK or SlimDX or SharpDX too) here is

			// slower! So it is disabled.

			return

				matrix1.M11 == matrix2.M11 &&

				matrix1.M22 == matrix2.M22 &&

				matrix1.M33 == matrix2.M33 &&

				matrix1.M44 == matrix2.M44 &&

				// Now check the rest if the diagonal matrices were the same.

				matrix1.M12 == matrix2.M12 &&

				matrix1.M13 == matrix2.M13 &&

				matrix1.M14 == matrix2.M14 &&

				matrix1.M21 == matrix2.M21 &&

				matrix1.M23 == matrix2.M23 &&

				matrix1.M24 == matrix2.M24 &&

				matrix1.M31 == matrix2.M31 &&

				matrix1.M32 == matrix2.M32 &&

				matrix1.M34 == matrix2.M34 &&

				matrix1.M41 == matrix2.M41 &&

				matrix1.M42 == matrix2.M42 &&

				matrix1.M43 == matrix2.M43;

		}

		#endregion



		#region Unequal (Static)

		/// <summary>

		/// Inequality check. Optimized for speed via ref parameters, which

		/// are much faster than copying 128 bytes (2 matrices each 64 bytes).

		/// </summary>

		public static bool Unequal(ref Matrix matrix1, ref Matrix matrix2)

		{

			// Check if everything is equal, if not return true (this is the unequal

			// check), might look a bit strange, but it is fast because of early out.

			if (matrix1.M11 == matrix2.M11 &&

			    matrix1.M12 == matrix2.M12 &&

			    matrix1.M13 == matrix2.M13 &&

			    matrix1.M14 == matrix2.M14 &&

			    matrix1.M21 == matrix2.M21 &&

			    matrix1.M22 == matrix2.M22 &&

			    matrix1.M23 == matrix2.M23 &&

			    matrix1.M24 == matrix2.M24 &&

			    matrix1.M31 == matrix2.M31 &&

			    matrix1.M32 == matrix2.M32 &&

			    matrix1.M33 == matrix2.M33 &&

			    matrix1.M34 == matrix2.M34 &&

			    matrix1.M41 == matrix2.M41 &&

			    matrix1.M42 == matrix2.M42 &&

			    matrix1.M43 == matrix2.M43)

			{

				return (matrix1.M44 != matrix2.M44);

			}



			return true;

		}

		#endregion



		#region TranslationNearlyEqual (Static)

		/// <summary>

		/// Nearly equals (only checks the Translation of the matrices), two

		/// matrices are considered nearly equal if their translation distance

		/// is below MathHelper.Epsilon. Use the MatrixNearlyEqual method to

		/// compare two complete matrices (much slower than this check).

		/// </summary>

		/// <param name="matrix1">Matrix 1 to compare</param>

		/// <param name="matrix2">Matrix 2 to compare</param>

		/// <returns>True if the translations of the matrices are almost the

		/// same, false otherwise.</returns>

		public static bool TranslationNearlyEqual(ref Matrix matrix1,

			ref Matrix matrix2)

		{

			// Check to see if the absolute difference of each value is smaller than

			// the Nearly Equal value.

			return (matrix1.Translation - matrix2.Translation).LengthSquared <

			       MathHelper.Epsilon * MathHelper.Epsilon;

		}

		#endregion



		#region MatrixNearlyEqual (Static)

		/// <summary>

		/// Matrix nearly equal helper method to compare two matrices that might

		/// not be 100% the same because of rounding errors, but they produce

		/// pretty much the same results.

		/// </summary>

		/// <param name="matrix1">Matrix 1 to compare</param>

		/// <param name="matrix2">Matrix 2 to compare</param>

		/// <returns>True if the matrices are almost the same</returns>

		public static bool MatrixNearlyEqual(ref Matrix mat1, ref Matrix mat2)

		{

			return

				mat1.M11.NearlyEqual(mat2.M11) &&

				mat1.M12.NearlyEqual(mat2.M12) &&

				mat1.M13.NearlyEqual(mat2.M13) &&

				mat1.M14.NearlyEqual(mat2.M14) &&

				mat1.M21.NearlyEqual(mat2.M21) &&

				mat1.M22.NearlyEqual(mat2.M22) &&

				mat1.M23.NearlyEqual(mat2.M23) &&

				mat1.M24.NearlyEqual(mat2.M24) &&

				mat1.M31.NearlyEqual(mat2.M31) &&

				mat1.M32.NearlyEqual(mat2.M32) &&

				mat1.M33.NearlyEqual(mat2.M33) &&

				mat1.M34.NearlyEqual(mat2.M34) &&

				mat1.M41.NearlyEqual(mat2.M41) &&

				mat1.M42.NearlyEqual(mat2.M42) &&

				mat1.M43.NearlyEqual(mat2.M43) &&

				mat1.M44.NearlyEqual(mat2.M44);

		}

		#endregion



		#region FromQuaternion (Static)

		/// <summary>

		/// From quaternion

		/// </summary>

		public static Matrix FromQuaternion(Quaternion quaternion)

		{

			//Check performance difference again:

			// Note: We can only be sure that this is faster on the Xbox, not used

			// often anyway.



			//Note: This could need some more comments, help, etc.

			float qxx = quaternion.X * quaternion.X;

			float qyy = quaternion.Y * quaternion.Y;

			float qzz = quaternion.Z * quaternion.Z;

			float qxy = quaternion.X * quaternion.Y;

			float qzw = quaternion.Z * quaternion.W;

			float qwx = quaternion.Z * quaternion.X;

			float qyw = quaternion.Y * quaternion.W;

			float qyz = quaternion.Y * quaternion.Z;

			float qxw = quaternion.X * quaternion.W;



			return new Matrix(

				1f - (2f * (qyy + qzz)), 2f * (qxy + qzw), 2f * (qwx - qyw), 0f,

				2f * (qxy - qzw), 1f - (2f * (qzz + qxx)), 2f * (qyz + qxw), 0f,

				2f * (qwx + qyw), 2f * (qyz - qxw), 1f - (2f * (qyy + qxx)), 0f,

				0f, 0f, 0f, 1f);

		}

		#endregion



		#region FromString (Static)

		/// <summary>

		/// Convert a string to a Matrix. The expected format is

		/// (M11, M12, M13, M14)

		/// (M21, M22, M23, M24)

		/// ...

		/// </summary>

		/// <param name="matrixString">The string containing the values in the

		/// correct format.</param>

		public static Matrix FromString(string matrixString)

		{

			// First remove the brackets

			matrixString = matrixString.

				Replace("(", "").

				Replace(")", "");

			// And split the string up into seperate values.

			string[] matrixStrings = matrixString.SplitAndTrim(',');

			//StringSplitOptions.RemoveEmptyEntries);



			// Then check if the length is 16 for 16 values and return the new matrix.

			// If the length is not 16 than return Matrix.Identity.

			if (matrixStrings.Length == 16)

			{

				return new Matrix(

					// First row

					matrixStrings[0].FromInvariantString(1.0f),

					matrixStrings[1].FromInvariantString(0.0f),

					matrixStrings[2].FromInvariantString(0.0f),

					matrixStrings[3].FromInvariantString(0.0f),

					// Second row

					matrixStrings[4].FromInvariantString(0.0f),

					matrixStrings[5].FromInvariantString(1.0f),

					matrixStrings[6].FromInvariantString(0.0f),

					matrixStrings[7].FromInvariantString(0.0f),

					// Third row

					matrixStrings[8].FromInvariantString(0.0f),

					matrixStrings[9].FromInvariantString(0.0f),

					matrixStrings[10].FromInvariantString(1.0f),

					matrixStrings[11].FromInvariantString(0.0f),

					// Forth row

					matrixStrings[12].FromInvariantString(0.0f),

					matrixStrings[13].FromInvariantString(0.0f),

					matrixStrings[14].FromInvariantString(0.0f),

					matrixStrings[15].FromInvariantString(1.0f));

			}



			Log.Warning("Unable to convert matrixString=" + matrixString +

			            " because it has not exactly 16 float values!");

			return Identity;

		}

		#endregion



		#region Identity (Static)

		/// <summary>

		/// Returns a new identity matrix (no rotation, translation or scaling)!

		/// Note: Not readonly to allow passing it by ref, but this should never

		/// be set!

		/// </summary>

		public static Matrix Identity = new Matrix(

			1, 0, 0, 0,

			0, 1, 0, 0,

			0, 0, 1, 0,

			0, 0, 0, 1);

		#endregion



		#region Framework Union Defines (Public)

		#endregion



		#region M11 (Public)

		/// <summary>

		/// 1st row - 1st column value of the Matrix.

		/// </summary>

		[FieldOffset(0)]

		public float M11;

		#endregion



		#region M12 (Public)

		/// <summary>

		/// 1st row - 2nd column value of the Matrix.

		/// </summary>

		[FieldOffset(4)]

		public float M12;

		#endregion



		#region M13 (Public)

		/// <summary>

		/// 1st row - 3rd column value of the Matrix.

		/// </summary>

		[FieldOffset(8)]

		public float M13;

		#endregion



		#region M14 (Public)

		/// <summary>

		/// 1st row - 4th column value of the Matrix.

		/// </summary>

		[FieldOffset(12)]

		public float M14;

		#endregion



		#region M21 (Public)

		/// <summary>

		/// 2nd row - 1st column value of the Matrix.

		/// </summary>

		[FieldOffset(16)]

		public float M21;

		#endregion



		#region M22 (Public)

		/// <summary>

		/// 2nd row - 2nd column value of the Matrix.

		/// </summary>

		[FieldOffset(20)]

		public float M22;

		#endregion



		#region M23 (Public)

		/// <summary>

		/// 2nd row - 3nd column value of the Matrix.

		/// </summary>

		[FieldOffset(24)]

		public float M23;

		#endregion



		#region M24 (Public)

		/// <summary>

		/// 2nd row - 4nd column value of the Matrix.

		/// </summary>

		[FieldOffset(28)]

		public float M24;

		#endregion



		#region M31 (Public)

		/// <summary>

		/// 3rd row - 1st column value of the Matrix.

		/// </summary>

		[FieldOffset(32)]

		public float M31;

		#endregion



		#region M32 (Public)

		/// <summary>

		/// 3rd row - 2nd column value of the Matrix.

		/// </summary>

		[FieldOffset(36)]

		public float M32;

		#endregion



		#region M33 (Public)

		/// <summary>

		/// 3rd row - 3rd column value of the Matrix.

		/// </summary>

		[FieldOffset(40)]

		public float M33;

		#endregion



		#region M34 (Public)

		/// <summary>

		/// 3rd row - 4th column value of the Matrix.

		/// </summary>

		[FieldOffset(44)]

		public float M34;

		#endregion



		#region M41 (Public)

		/// <summary>

		/// 4th row - 1st column value of the Matrix.

		/// </summary>

		[FieldOffset(48)]

		public float M41;

		#endregion



		#region M42 (Public)

		/// <summary>

		/// 4th row - 2nd column value of the Matrix.

		/// </summary>

		[FieldOffset(52)]

		public float M42;

		#endregion



		#region M43 (Public)

		/// <summary>

		/// 4th row - 3rd column value of the Matrix.

		/// </summary>

		[FieldOffset(56)]

		public float M43;

		#endregion



		#region M44 (Public)

		/// <summary>

		/// 4th row - 4th column value of the Matrix.

		/// </summary>

		[FieldOffset(60)]

		public float M44;

		#endregion



		#region Translation (Public)

		/// <summary>

		/// The translation amount representing by the matrix.

		/// This vector shares the same data as M41, M42, M43, see above!

		/// </summary>

		[FieldOffset(48)]

		public Vector Translation;

		#endregion



		#region Scaling (Public)

		/// <summary>

		/// The scaling amount representing by the matrix.

		/// </summary>

		public Vector Scaling

		{

			get

			{

				return new Vector(M11, M22, M33);

			}

			set

			{

				M11 = value.X;

				M22 = value.Y;

				M33 = value.Z;

			}

		}

		#endregion



		#region Right (Public)

		/// <summary>

		/// Right vector of the matrix (M11, M12, M13)

		/// </summary>

		[FieldOffset(0)]

		public Vector Right;

		#endregion



		#region Up (Public)

		/// <summary>

		/// Up vector of the matrix (M21, M22, M23)

		/// </summary>

		[FieldOffset(16)]

		public Vector Up;

		#endregion



		#region Front (Public)

		/// <summary>

		/// Front vector of the matrix (M31, M32, M33)

		/// </summary>

		[FieldOffset(32)]

		public Vector Front;

		#endregion



		#region Determinant (Public)

		/// <summary>

		/// Returns the determinant of the matrix. Uses optimized code paths

		/// because XNAs implementation is twice as fast as ours, OpenTK is

		/// slightly faster (see MatrixPerformance class).

		/// </summary>

		public float Determinant

		{

			get

			{

				//Check performance difference again:



				// Ugly hacky code that is actually lots faster than our clean

				// implementation below (see link on how it works basically).

				float m44x33m43x34 = (M44 * M33) - (M43 * M34);

				float m32x44m42x34 = (M32 * M44) - (M42 * M34);

				float m32x43m42x33 = (M32 * M43) - (M42 * M33);

				float m31x44m41x34 = (M31 * M44) - (M41 * M34);

				float m31x43m41x33 = (M31 * M43) - (M41 * M33);

				float m31x42m41x32 = (M31 * M42) - (M41 * M32);

				return (((((((M22 * m44x33m43x34) - (M23 * m32x44m42x34)) +

				            (M24 * m32x43m42x33)) * M11) - ((((M21 * m44x33m43x34) -

				                                              (M23 * m31x44m41x34)) +

				                                             (M24 * m31x43m41x33)) * M12)) +

				         ((((M21 * m32x44m42x34) - (M22 * m31x44m41x34)) +

				           (M24 * m31x42m41x32)) * M13)) - ((((M21 * m32x43m42x33) -

				                                              (M22 * m31x43m41x33)) +

				                                             (M23 * m31x42m41x32)) * M14));

			}

		}

		#endregion



		#region Inverse (Public)

		/// <summary>

		/// Returns the inverse of the current matrix, which can be slightly

		/// faster than using the static Invert method below depending on the

		/// platform (e.g. with XNA), with Delta matrix code Invert() is a little

		/// faster however, so just use whatever method fits best for you.

		/// Both methods are the slowest matrix methods however, use with care!

		/// </summary>

		public Matrix Inverse

		{

			get

			{

				//Check performance difference again:

				// Ugly hacky code that is actually lots faster than our clean

				// implementation below (see link on how it works basically).

				// This is basically getting the Determinant and then also calculating

				// the inverse matrix in one big ugly swoop.

				float m33x44m34x43 = (M33 * M44) - (M34 * M43);

				float m32x44m34x42 = (M32 * M44) - (M34 * M42);

				float m32x43m33x42 = (M32 * M43) - (M33 * M42);

				float m31x44m34x41 = (M31 * M44) - (M34 * M41);

				float m31x43m33x41 = (M31 * M43) - (M33 * M41);

				float m31x42m32x41 = (M31 * M42) - (M32 * M41);

				float m22m23m24 = ((M22 * m33x44m34x43) - (M23 * m32x44m34x42)) +

				                  (M24 * m32x43m33x42);

				float m21m23m24 = -(((M21 * m33x44m34x43) - (M23 * m31x44m34x41)) +

				                    (M24 * m31x43m33x41));

				float m21m22m24 = ((M21 * m32x44m34x42) - (M22 * m31x44m34x41)) +

				                  (M24 * m31x42m32x41);

				float m21m22m23 = -(((M21 * m32x43m33x42) - (M22 * m31x43m33x41)) +

				                    (M23 * m31x42m32x41));

				// Inverse(A) = 1/det(A) * B

				float inverseDet = 1f / ((((M11 * m22m23m24) + (M12 * m21m23m24)) +

				                          (M13 * m21m22m24)) + (M14 * m21m22m23));

				float m23x44m24x43 = (M23 * M44) - (M24 * M43);

				float m22x44m24x42 = (M22 * M44) - (M24 * M42);

				float m22x43m23x42 = (M22 * M43) - (M23 * M42);

				float m21x44m24x41 = (M21 * M44) - (M24 * M41);

				float m21x43m23x41 = (M21 * M43) - (M23 * M41);

				float m21x42m22x41 = (M21 * M42) - (M22 * M41);

				float m23x34m24x33 = (M23 * M34) - (M24 * M33);

				float m22x34m24x32 = (M22 * M34) - (M24 * M32);

				float m22x33m23x32 = (M22 * M33) - (M23 * M32);

				float m21x34m23x31 = (M21 * M34) - (M24 * M31);

				float m21x33m23x31 = (M21 * M33) - (M23 * M31);

				float m21x32m22x31 = (M21 * M32) - (M22 * M31);

				return new Matrix(

					m22m23m24 * inverseDet,

					-(((M12 * m33x44m34x43) - (M13 * m32x44m34x42)) +

					  (M14 * m32x43m33x42)) * inverseDet,

					(((M12 * m23x44m24x43) - (M13 * m22x44m24x42)) +

					 (M14 * m22x43m23x42)) * inverseDet,

					-(((M12 * m23x34m24x33) - (M13 * m22x34m24x32)) +

					  (M14 * m22x33m23x32)) * inverseDet,

					m21m23m24 * inverseDet,

					(((M11 * m33x44m34x43) - (M13 * m31x44m34x41)) +

					 (M14 * m31x43m33x41)) * inverseDet,

					-(((M11 * m23x44m24x43) - (M13 * m21x44m24x41)) +

					  (M14 * m21x43m23x41)) * inverseDet,

					(((M11 * m23x34m24x33) - (M13 * m21x34m23x31)) +

					 (M14 * m21x33m23x31)) * inverseDet,

					m21m22m24 * inverseDet,

					-(((M11 * m32x44m34x42) - (M12 * m31x44m34x41)) +

					  (M14 * m31x42m32x41)) * inverseDet,

					(((M11 * m22x44m24x42) - (M12 * m21x44m24x41)) +

					 (M14 * m21x42m22x41)) * inverseDet,

					-(((M11 * m22x34m24x32) - (M12 * m21x34m23x31)) +

					  (M14 * m21x32m22x31)) * inverseDet,

					m21m22m23 * inverseDet,

					(((M11 * m32x43m33x42) - (M12 * m31x43m33x41)) +

					 (M13 * m31x42m32x41)) * inverseDet,

					-(((M11 * m22x43m23x42) - (M12 * m21x43m23x41)) +

					  (M13 * m21x42m22x41)) * inverseDet,

					(((M11 * m22x33m23x32) - (M12 * m21x33m23x31)) +

					 (M13 * m21x32m22x31)) * inverseDet);

			}

		}

		#endregion



		#region Constructors

		/// <summary>

		/// Create matrix by passing all the necessary values

		/// </summary>

		/// <param name="setM11">Set M11 value (first row)</param>

		/// <param name="setM12">Set M12 value (first row)</param>

		/// <param name="setM13">Set M13 value (first row)</param>

		/// <param name="setM14">Set M14 value (first row)</param>

		/// <param name="setM21">Set M21 value (second row)</param>

		/// <param name="setM22">Set M22 value (second row)</param>

		/// <param name="setM23">Set M23 value (second row)</param>

		/// <param name="setM24">Set M24 value (second row)</param>

		/// <param name="setM31">Set M31 value (third row)</param>

		/// <param name="setM32">Set M32 value (third row)</param>

		/// <param name="setM33">Set M33 value (third row)</param>

		/// <param name="setM34">Set M34 value (third row)</param>

		/// <param name="setM41">Set M41 value (forth row)</param>

		/// <param name="setM42">Set M42 value (forth row)</param>

		/// <param name="setM43">Set M43 value (forth row)</param>

		/// <param name="setM44">Set M44 value (forth row)</param>

		public Matrix(float setM11, float setM12, float setM13, float setM14,

			float setM21, float setM22, float setM23, float setM24,

			float setM31, float setM32, float setM33, float setM34,

			float setM41, float setM42, float setM43, float setM44)

			: this()

		{

			M11 = setM11;

			M12 = setM12;

			M13 = setM13;

			M14 = setM14;

			M21 = setM21;

			M22 = setM22;

			M23 = setM23;

			M24 = setM24;

			M31 = setM31;

			M32 = setM32;

			M33 = setM33;

			M34 = setM34;

			M41 = setM41;

			M42 = setM42;

			M43 = setM43;

			M44 = setM44;

		}



		/// <summary>

		/// Create matrix by passing all 16 matrix values as array.

		/// m11, m12, m13, m14, 

		/// m21, m22, m23, m24,

		/// m31, m32, m33, m34, 

		/// m41, m42, m43, m44, 

		/// </summary>

		/// <param name="setMatrix">Set Matrix</param>

		public Matrix(float[] setMatrix)

			: this()

		{

			int length = setMatrix.Length;

			if (length < 16)

			{

				Log.Warning("The matrix has too little values.\n" +

				            "Length: " + length);

			}

			if (length > 16)

			{

				Log.Warning("The matrix has too many values.\n" +

				            "Length: " + length);

			}



			// Copy all values.

			M11 = setMatrix[0];

			M12 = setMatrix[1];

			M13 = setMatrix[2];

			M14 = setMatrix[3];

			M21 = setMatrix[4];

			M22 = setMatrix[5];

			M23 = setMatrix[6];

			M24 = setMatrix[7];

			M31 = setMatrix[8];

			M32 = setMatrix[9];

			M33 = setMatrix[10];

			M34 = setMatrix[11];

			M41 = setMatrix[12];

			M42 = setMatrix[13];

			M43 = setMatrix[14];

			M44 = setMatrix[15];

		}



		/// <summary>

		/// Create matrix by passing in 3 vectors for M11-M33, rest stays default.

		/// </summary>

		/// <param name="firstRowVector">Vector for the first row (M11-M13)

		/// </param>

		/// <param name="secondRowVector">Vector for the second row (M21-M23)

		/// </param>

		/// <param name="thirdRowVector">Vector for the third row (M31-M33)

		/// </param>

		public Matrix(Vector firstRowVector, Vector secondRowVector,

			Vector thirdRowVector)

			: this()

		{

			// Copy all values and fill in the blanks

			M11 = firstRowVector.X;

			M12 = firstRowVector.Y;

			M13 = firstRowVector.Z;

			M14 = 0.0f;

			M21 = secondRowVector.X;

			M22 = secondRowVector.Y;

			M23 = secondRowVector.Z;

			M24 = 0.0f;

			M31 = thirdRowVector.X;

			M32 = thirdRowVector.Y;

			M33 = thirdRowVector.Z;

			M34 = 0.0f;

			M41 = 0.0f;

			M42 = 0.0f;

			M43 = 0.0f;

			M44 = 1.0f;

		}



		/// <summary>

		/// Create matrix with help of BinaryReader, we will read 16 floats from

		/// the stream that is linked up with this BinaryReader.

		/// </summary>

		/// <param name="reader">Reader</param>

		public Matrix(BinaryReader reader)

			: this()

		{

			Load(reader);

		}

		#endregion



		#region IEquatable<Matrix> Members

		/// <summary>

		/// Equals, quickly checks if another matrix has the exact same values.

		/// This methods needs to copy 64 bytes around (1 matrix of 64 bytes).

		/// Try to use the ref version instead.

		/// </summary>

		/// <param name="other">Other Matrix to compare to</param>

		/// <returns>True if both matrices are equal, false otherwise.</returns>

		public bool Equals(Matrix other)

		{

			// First check the diagonal values, those mostly differ.

			return

				M11 == other.M11 &&

				M22 == other.M22 &&

				M33 == other.M33 &&

				M44 == other.M44 &&

				// Now check the rest if the diagonal values were the same.

				M12 == other.M12 &&

				M13 == other.M13 &&

				M14 == other.M14 &&

				M21 == other.M21 &&

				M23 == other.M23 &&

				M24 == other.M24 &&

				M31 == other.M31 &&

				M32 == other.M32 &&

				M34 == other.M34 &&

				M41 == other.M41 &&

				M42 == other.M42 &&

				M43 == other.M43;

		}

		#endregion



		#region ISaveLoadBinary Members

		/// <summary>

		/// Load the matrix values from a stream (64 bytes, 16 floats).

		/// </summary>

		/// <param name="reader">The stream that will be used.</param>

		public void Load(BinaryReader reader)

		{

			// Note XNA and OpenTK matrices share the same data, no need to set them

			M11 = reader.ReadSingle();

			M12 = reader.ReadSingle();

			M13 = reader.ReadSingle();

			M14 = reader.ReadSingle();



			M21 = reader.ReadSingle();

			M22 = reader.ReadSingle();

			M23 = reader.ReadSingle();

			M24 = reader.ReadSingle();



			M31 = reader.ReadSingle();

			M32 = reader.ReadSingle();

			M33 = reader.ReadSingle();

			M34 = reader.ReadSingle();



			M41 = reader.ReadSingle();

			M42 = reader.ReadSingle();

			M43 = reader.ReadSingle();

			M44 = reader.ReadSingle();

		}



		/// <summary>

		/// Saves the matrix to a stream (64 bytes, 16 floats).

		/// </summary>

		/// <param name="writer">The stream that will be used.</param>

		public void Save(BinaryWriter writer)

		{

			writer.Write(M11);

			writer.Write(M12);

			writer.Write(M13);

			writer.Write(M14);



			writer.Write(M21);

			writer.Write(M22);

			writer.Write(M23);

			writer.Write(M24);



			writer.Write(M31);

			writer.Write(M32);

			writer.Write(M33);

			writer.Write(M34);



			writer.Write(M41);

			writer.Write(M42);

			writer.Write(M43);

			writer.Write(M44);

		}

		#endregion



		#region op_Multiply (Operator)

		/// <summary>

		/// Multiply operator. Warning: This operator is slower than using ref

		/// version, it needs to copy 192 bytes around (3 matrices, each 64 bytes).

		/// </summary>

		/// <param name="matrix1">Matrix 1 to multiply</param>

		/// <param name="matrix2">Matrix 2 to multiply</param>

		public static Matrix operator *(Matrix matrix1, Matrix matrix2)

		{

			// Other discussions about XNA Multiply performance

			// http://forums.xna.com/forums/t/7061.aspx?PageIndex=2

			// Machine: Core 2 Duo E6600 clocked at 3.0GHz

			// Managed XNA - Matrix.Multiply(ref,ref,out): 557ms

			// Managed w/ P/Invoke to SSE2 function:       341ms

			// Native C++ - Full Compiler Optimizations:   300ms

			// Native C++ - hand-optimized SSE2:           138ms

			// Note: We got like 97ms with this C# code, but a faster PC ^^



			// Try to do some early out code in case matrix2 is an identity matrix

			if (matrix2.M11 == 1.0f &&

			    matrix2.M22 == 1.0f &&

			    matrix2.M33 == 1.0f &&

			    matrix2.M44 == 1.0f &&

			    // Ok, just to be sure, we need to check if all the other values

			    // are really 0, only then this is an identity matrix

			    matrix2.M21 == 0.0f &&

			    matrix2.M23 == 0.0f &&

			    matrix2.M24 == 0.0f &&

			    matrix2.M31 == 0.0f &&

			    matrix2.M32 == 0.0f &&

			    matrix2.M34 == 0.0f &&

			    matrix2.M41 == 0.0f &&

			    matrix2.M42 == 0.0f &&

			    matrix2.M43 == 0.0f)

			{

				return matrix1;

			}



			return new Matrix(

				// M11

				(matrix1.M11 * matrix2.M11) + (matrix1.M12 * matrix2.M21) +

				(matrix1.M13 * matrix2.M31) + (matrix1.M14 * matrix2.M41),

				// M12

				(matrix1.M11 * matrix2.M12) + (matrix1.M12 * matrix2.M22) +

				(matrix1.M13 * matrix2.M32) + (matrix1.M14 * matrix2.M42),

				// M13

				(matrix1.M11 * matrix2.M13) + (matrix1.M12 * matrix2.M23) +

				(matrix1.M13 * matrix2.M33) + (matrix1.M14 * matrix2.M43),

				// M14

				(matrix1.M11 * matrix2.M14) + (matrix1.M12 * matrix2.M24) +

				(matrix1.M13 * matrix2.M34) + (matrix1.M14 * matrix2.M44),

				// M21

				(matrix1.M21 * matrix2.M11) + (matrix1.M22 * matrix2.M21) +

				(matrix1.M23 * matrix2.M31) + (matrix1.M24 * matrix2.M41),

				// M22

				(matrix1.M21 * matrix2.M12) + (matrix1.M22 * matrix2.M22) +

				(matrix1.M23 * matrix2.M32) + (matrix1.M24 * matrix2.M42),

				// M23

				(matrix1.M21 * matrix2.M13) + (matrix1.M22 * matrix2.M23) +

				(matrix1.M23 * matrix2.M33) + (matrix1.M24 * matrix2.M43),

				// M24

				(matrix1.M21 * matrix2.M14) + (matrix1.M22 * matrix2.M24) +

				(matrix1.M23 * matrix2.M34) + (matrix1.M24 * matrix2.M44),

				// M31

				(matrix1.M31 * matrix2.M11) + (matrix1.M32 * matrix2.M21) +

				(matrix1.M33 * matrix2.M31) + (matrix1.M34 * matrix2.M41),

				// M32

				(matrix1.M31 * matrix2.M12) + (matrix1.M32 * matrix2.M22) +

				(matrix1.M33 * matrix2.M32) + (matrix1.M34 * matrix2.M42),

				// M33

				(matrix1.M31 * matrix2.M13) + (matrix1.M32 * matrix2.M23) +

				(matrix1.M33 * matrix2.M33) + (matrix1.M34 * matrix2.M43),

				// M34

				(matrix1.M31 * matrix2.M14) + (matrix1.M32 * matrix2.M24) +

				(matrix1.M33 * matrix2.M34) + (matrix1.M34 * matrix2.M44),

				// M41

				(matrix1.M41 * matrix2.M11) + (matrix1.M42 * matrix2.M21) +

				(matrix1.M43 * matrix2.M31) + (matrix1.M44 * matrix2.M41),

				// M42

				(matrix1.M41 * matrix2.M12) + (matrix1.M42 * matrix2.M22) +

				(matrix1.M43 * matrix2.M32) + (matrix1.M44 * matrix2.M42),

				// M43

				(matrix1.M41 * matrix2.M13) + (matrix1.M42 * matrix2.M23) +

				(matrix1.M43 * matrix2.M33) + (matrix1.M44 * matrix2.M43),

				// M44

				(matrix1.M41 * matrix2.M14) + (matrix1.M42 * matrix2.M24) +

				(matrix1.M43 * matrix2.M34) + (matrix1.M44 * matrix2.M44));

		}



		/// <summary>

		/// Operator to multiply matrix with a vector. This is very useful to

		/// transform vectors with the given matrix. Warning: This is slower than

		/// the ref version, this operator needs to copy 88 bytes around (1 matrix,

		/// 64 bytes and 2 vectors, each 12 bytes)

		/// </summary>

		/// <param name="matrix">Matrix to multiply</param>

		/// <param name="vector">Vector to multiply with</param>

		/// <returns>Resulting vector from the multiplication.</returns>

		public static Vector operator *(Matrix matrix, Vector vector)

		{

			// Note: Seems like calling XNA (OpenTK or SlimDX too?) is here slower!

			//Delta: 209ms (10 million calls)

			//XNA: 353ms (10 million calls)

			//OpenTK: 245ms (10 million calls)

			//SlimDX: 265ms (10 million calls)

			// Just use the normal vector * matrix formula (see Vector class)

			return new Vector(

				(vector.X * matrix.M11) + (vector.Y * matrix.M21) +

				(vector.Z * matrix.M31) + matrix.M41,

				(vector.X * matrix.M12) + (vector.Y * matrix.M22) +

				(vector.Z * matrix.M32) + matrix.M42,

				(vector.X * matrix.M13) + (vector.Y * matrix.M23) +

				(vector.Z * matrix.M33) + matrix.M43);

		}



		/// <summary>

		/// Operator to multiply matrix with a scale factor. Try to use the ref

		/// version instead, this needs to copy 128 bytes around (2 matrices, each

		/// 64 bytes).

		/// </summary>

		/// <param name="matrix">Matrix to multiply</param>

		/// <param name="scaleFactor">Scale factor to multiply with</param>

		/// <returns>

		/// New matrix with all values multiplied by scaleFactor.

		/// </returns>

		public static Matrix operator *(Matrix matrix, float scaleFactor)

		{

			//Note: XNA, OpenTK, and SlimDX path