C#: Camera Matrix

I always wanted to know how a 3D engine works. I've used DirectX several times before, but I never quite understood what's going on behind the scene.
So, I decided to write a small naive 3D engine of my own, just to understand how it works.
I studied about the subject, and wrote this engine, which is basically a camera matrix: (no textures or anything)
http://www.megaupload.com/?d=6YON627Z
The only external graphics function I used is the Graphics.DrawLine() function.

What is a Camera Matrix anyway?
It's something that can be used to display 3D objects on our 2D screen.
Think about a real camera like in a 2D movie (not Avatar) - you place it somewhere in the room (translation), it looks towards the scene (rotation), it can zoom in and out (scaling), and it converts the 3D scene to 2D (projection).
From now on we say "translation" instead of "position". Mathematicians are to blame for that.

It's easy to do operations on vectors, why complicate things with Matrices?
One of the main reasons is that matrices can be easily concatenated.
From DirectX MSDN:

One advantage of using matrices is that you can combine the effects of two or more matrices by multiplying them. This means that, to rotate a model and then translate it to some location, you do not need to apply two matrices. Instead, you multiply the rotation and translation matrices to produce a composite matrix that contains all of their effects.

Let me put it this way:
"Using" a 4x4 matrix requires 16 multiplications and 9 additions.
So rotating, which is "using" a rotation matrix, requires 16 multiplications and 9 additions.
And scaling also requires 16 multiplications and 9 additions.
If we do rotating and then scaling it would take 32 multiplications and 18 additions. Translation would takes us up to 48 multiplications and 27 additions, and so on...
BUT, if we concatenate the rotation and scaling matrices, we will get a third matrix that does both operations in half the time (16 multiplications and 9 additions).
In fact, we could concatenate any number of matrices, so the resulting matrix will do all the operations we want in only 16 multiplications and 9 additions.
You might have noticed that we still need to do the concatenation operation which takes a long time, but it's worth it when you design a 3D game with thousands and thousands of objects to scale rotate translate etc. but only a few matrices to concatenate.

Why a Camera Matrix is not 3x3?
In short, because then we can't use it for translation (moving the vectors).
For more details read this fantastic article:
http://www.geometer.org/mathcircles/cghomogen.pdf

C#: Matrix and Vector Representation

In C# there are built in representations of 2D points and 2x2 matrices.
This is good enough for image processing, but I needed more than that for 3D graphics.

I wrote classes for:

• 3D vector
• 3x3 matrix
• 4x4 matrix
• Any size matrix

All the first three objects can be represented by "Any size matrix", but it is much more efficient to use the specific class you need.

Code:

class Vector3
{
    public static Vector3 Zero = NewZero();
    public static Vector3 One = NewOne();

    public float x;
    public float y;
    public float z;

    public Vector3()
    {
        x = 0.0f;
        y = 0.0f;
        z = 0.0f;
    }

    public Vector3(float x, float y, float z)
    {
        this.x = x;
        this.y = y;
        this.z = z;
    }

    public Vector3(float xyz)
    {
        this.x = xyz;
        this.y = xyz;
        this.z = xyz;
    }

    public Vector3(Vector3 v)
    {
        this.x = v.x;
        this.y = v.y;
        this.z = v.z;
    }

    public static Vector3 NewZero()
    {
        return new Vector3(0.0f);
    }

    public static Vector3 NewOne()
    {
        return new Vector3(1.0f);
    }

    public float DotProduct(Vector3 other)
    {
        return x * other.x + y * other.y + z * other.z;
    }

    public static Vector3 operator +(Vector3 v1, Vector3 v2)
    {
        return new Vector3(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z);
    }

    public static Vector3 operator -(Vector3 v1, Vector3 v2)
    {
        return new Vector3(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z);
    }

    public static Vector3 operator -(Vector3 v)
    {
        return new Vector3(-v.x, -v.y, -v.z);
    }

    public static Vector3 operator *(Vector3 v, float scalar)
    {
        return new Vector3(v.x * scalar, v.y * scalar, v.z * scalar);
    }

    public static Vector3 operator /(Vector3 v, float scalar)
    {
        return new Vector3(v.x / scalar, v.y / scalar, v.z / scalar);
    }

    public static bool operator ==(Vector3 v1, Vector3 v2)
    {
        return v1.x == v2.x && v1.y == v2.y && v1.z == v2.z;
    }

    public static bool operator !=(Vector3 v1, Vector3 v2)
    {
        return v1.x != v2.x || v1.y != v2.y || v1.z != v2.z;
    }

    public static Vector3 CrossProduct(Vector3 a, Vector3 b)
    {
        return new Vector3(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x);
    }

    public Vector3 Add(Vector3 v)
    {
        x += v.x;
        y += v.y;
        z += v.z;
        return this;
    }

    public float DistanceTo(Vector3 v)
    {
        float dx = this.x - v.x;
        float dy = this.y - v.y;
        float dz = this.z - v.z;
        return (float)Math.Sqrt(dx * dx + dy * dy + dz * dz);
    }

    public float Size()
    {
        return DistanceTo(Vector3.Zero);
    }

    public Vector3 Normalize()
    {
        float size = Size();
        this.x /= size;
        this.y /= size;
        this.z /= size;
        return this;
    }

    public Vector3 Clone()
    {
        return new Vector3(this);
    }

    public override string ToString()
    {
        return "(" + x + ", " + y + ", " + z + ")";
    }
}

class Matrix
{
    public float[,] matrix;
    public int rows;
    public int cols;

    public Matrix(int rows, int cols)
    {
        this.matrix = new float[rows, cols];
        this.rows = rows;
        this.cols = cols;
    }

    public Matrix(float[,] matrix)
    {
        this.matrix = matrix;
        this.rows = matrix.GetLength(0);
        this.cols = matrix.GetLength(1);
    }

    protected static float[,] Multiply(Matrix matrix, float scalar)
    {
        int rows = matrix.rows;
        int cols = matrix.cols;
        float[,] m1 = matrix.matrix;
        float[,] m2 = new float[rows, cols];
        for (int i = 0; i < rows; ++i)
        {
            for (int j = 0; j < cols; ++j)
            {
                m2[i, j] = m1[i, j] * scalar;
            }
        }
        return m2;
    }

    protected static float[,] Multiply(Matrix matrix1, Matrix matrix2)
    {
        int m1rows = matrix1.rows;
        int m1cols = matrix1.cols;
        int m2rows = matrix2.rows;
        int m2cols = matrix2.cols;
        if (m1cols != m2rows)
        {
            throw new ArgumentException();
        }
        float[,] m1 = matrix1.matrix;
        float[,] m2 = matrix2.matrix;
        float[,] m3 = new float[m1rows, m2cols];
        for (int i = 0; i < m1rows; ++i)
        {
            for (int j = 0; j < m2cols; ++j)
            {
                float sum = 0;
                for (int it = 0; it < m1cols; ++it)
                {
                    sum += m1[i, it] * m2[it, j];
                }
                m3[i, j] = sum;
            }
        }
        return m3;
    }

    public static Matrix operator *(Matrix m, float scalar)
    {
        return new Matrix(Multiply(m, scalar));
    }

    public static Matrix operator *(Matrix m1, Matrix m2)
    {
        return new Matrix(Multiply(m1, m2));
    }

    public override string ToString()
    {
        string res = "";
        for (int i = 0; i < rows; ++i)
        {
            if (i > 0)
            {
                res += "|";
            }
            for (int j = 0; j < cols; ++j)
            {
                if (j > 0)
                {
                    res += ",";
                }
                res += matrix[i, j];
            }
        }
        return "(" + res + ")";
    }
}

class Matrix3 : Matrix
{
    public Matrix3()
        : base(3, 3)
    {
    }

    public Matrix3(float[,] matrix)
        : base(matrix)
    {
        if (rows != 3 || cols != 3)
        {
            throw new ArgumentException();
        }
    }

    public static Matrix3 I()
    {
        return new Matrix3(new float[,] { 
        { 1.0f, 0.0f, 0.0f }, 
        { 0.0f, 1.0f, 0.0f }, 
        { 0.0f, 0.0f, 1.0f } });
    }

    public static Vector3 operator *(Matrix3 matrix3, Vector3 v)
    {
        float[,] m = matrix3.matrix;
        return new Vector3(
            m[0, 0] * v.x + m[0, 1] * v.y + m[0, 2] * v.z,
            m[1, 0] * v.x + m[1, 1] * v.y + m[1, 2] * v.z,
            m[2, 0] * v.x + m[2, 1] * v.y + m[2, 2] * v.z);
    }

    public static Matrix3 operator *(Matrix3 mat1, Matrix3 mat2)
    {
        float[,] m1 = mat1.matrix;
        float[,] m2 = mat2.matrix;
        float[,] m3 = new float[3, 3];
        m3[0, 0] = m1[0, 0] * m2[0, 0] + m1[0, 1] * m2[1, 0] + m1[0, 2] * m2[2, 0];
        m3[0, 1] = m1[0, 0] * m2[0, 1] + m1[0, 1] * m2[1, 1] + m1[0, 2] * m2[2, 1];
        m3[0, 2] = m1[0, 0] * m2[0, 2] + m1[0, 1] * m2[1, 2] + m1[0, 2] * m2[2, 2];
        m3[1, 0] = m1[1, 0] * m2[0, 0] + m1[1, 1] * m2[1, 0] + m1[1, 2] * m2[2, 0];
        m3[1, 1] = m1[1, 0] * m2[0, 1] + m1[1, 1] * m2[1, 1] + m1[1, 2] * m2[2, 1];
        m3[1, 2] = m1[1, 0] * m2[0, 2] + m1[1, 1] * m2[1, 2] + m1[1, 2] * m2[2, 2];
        m3[2, 0] = m1[2, 0] * m2[0, 0] + m1[2, 1] * m2[1, 0] + m1[2, 2] * m2[2, 0];
        m3[2, 1] = m1[2, 0] * m2[0, 1] + m1[2, 1] * m2[1, 1] + m1[2, 2] * m2[2, 1];
        m3[2, 2] = m1[2, 0] * m2[0, 2] + m1[2, 1] * m2[1, 2] + m1[2, 2] * m2[2, 2];
        return new Matrix3(m3);
    }

    public static Matrix3 operator *(Matrix3 m, float scalar)
    {
        return new Matrix3(Multiply(m, scalar));
    }
}

class Matrix4 : Matrix
{
    public static Matrix4 I = NewI();

    public Matrix4()
        : base(4, 4)
    {
    }

    public Matrix4(float[,] matrix)
        : base(matrix)
    {
        if (rows != 4 || cols != 4)
        {
            throw new ArgumentException();
        }
    }

    public static Matrix4 NewI()
    {
        return new Matrix4(new float[,] { 
        { 1.0f, 0.0f, 0.0f, 0.0f }, 
        { 0.0f, 1.0f, 0.0f, 0.0f }, 
        { 0.0f, 0.0f, 1.0f, 0.0f },
        { 0.0f, 0.0f, 0.0f, 1.0f } });
    }

    public static Vector3 operator *(Matrix4 matrix4, Vector3 v)
    {
        float[,] m = matrix4.matrix;
        float w = m[3, 0] * v.x + m[3, 1] * v.y + m[3, 2] * v.z + m[3, 3];
        return new Vector3(
            (m[0, 0] * v.x + m[0, 1] * v.y + m[0, 2] * v.z + m[0, 3]) / w,
            (m[1, 0] * v.x + m[1, 1] * v.y + m[1, 2] * v.z + m[1, 3]) / w,
            (m[2, 0] * v.x + m[2, 1] * v.y + m[2, 2] * v.z + m[2, 3]) / w
            );
    }

    public static Matrix4 operator *(Matrix4 mat1, Matrix4 mat2)
    {
        float[,] m1 = mat1.matrix;
        float[,] m2 = mat2.matrix;
        float[,] m3 = new float[4, 4];
        m3[0, 0] = m1[0, 0] * m2[0, 0] + m1[0, 1] * m2[1, 0] + m1[0, 2] * m2[2, 0] + m1[0, 3] * m2[3, 0];
        m3[0, 1] = m1[0, 0] * m2[0, 1] + m1[0, 1] * m2[1, 1] + m1[0, 2] * m2[2, 1] + m1[0, 3] * m2[3, 1];
        m3[0, 2] = m1[0, 0] * m2[0, 2] + m1[0, 1] * m2[1, 2] + m1[0, 2] * m2[2, 2] + m1[0, 3] * m2[3, 2];
        m3[0, 3] = m1[0, 0] * m2[0, 3] + m1[0, 1] * m2[1, 3] + m1[0, 2] * m2[2, 3] + m1[0, 3] * m2[3, 3];
        m3[1, 0] = m1[1, 0] * m2[0, 0] + m1[1, 1] * m2[1, 0] + m1[1, 2] * m2[2, 0] + m1[1, 3] * m2[3, 0];
        m3[1, 1] = m1[1, 0] * m2[0, 1] + m1[1, 1] * m2[1, 1] + m1[1, 2] * m2[2, 1] + m1[1, 3] * m2[3, 1];
        m3[1, 2] = m1[1, 0] * m2[0, 2] + m1[1, 1] * m2[1, 2] + m1[1, 2] * m2[2, 2] + m1[1, 3] * m2[3, 2];
        m3[1, 3] = m1[1, 0] * m2[0, 3] + m1[1, 1] * m2[1, 3] + m1[1, 2] * m2[2, 3] + m1[1, 3] * m2[3, 3];
        m3[2, 0] = m1[2, 0] * m2[0, 0] + m1[2, 1] * m2[1, 0] + m1[2, 2] * m2[2, 0] + m1[2, 3] * m2[3, 0];
        m3[2, 1] = m1[2, 0] * m2[0, 1] + m1[2, 1] * m2[1, 1] + m1[2, 2] * m2[2, 1] + m1[2, 3] * m2[3, 1];
        m3[2, 2] = m1[2, 0] * m2[0, 2] + m1[2, 1] * m2[1, 2] + m1[2, 2] * m2[2, 2] + m1[2, 3] * m2[3, 2];
        m3[2, 3] = m1[2, 0] * m2[0, 3] + m1[2, 1] * m2[1, 3] + m1[2, 2] * m2[2, 3] + m1[2, 3] * m2[3, 3];
        m3[3, 0] = m1[3, 0] * m2[0, 0] + m1[3, 1] * m2[1, 0] + m1[3, 2] * m2[2, 0] + m1[3, 3] * m2[3, 0];
        m3[3, 1] = m1[3, 0] * m2[0, 1] + m1[3, 1] * m2[1, 1] + m1[3, 2] * m2[2, 1] + m1[3, 3] * m2[3, 1];
        m3[3, 2] = m1[3, 0] * m2[0, 2] + m1[3, 1] * m2[1, 2] + m1[3, 2] * m2[2, 2] + m1[3, 3] * m2[3, 2];
        m3[3, 3] = m1[3, 0] * m2[0, 3] + m1[3, 1] * m2[1, 3] + m1[3, 2] * m2[2, 3] + m1[3, 3] * m2[3, 3];
        return new Matrix4(m3);
    }

    public static Matrix4 operator *(Matrix4 m, float scalar)
    {
        return new Matrix4(Multiply(m, scalar));
    }
}

C#: Rotation Matrix

Rotation matrices are used in 3D graphics to rotate vectors.
I use homogeneous coordinates, so the matrices are 4x4. If you want 3x3, just remove the last column and last row.
I wrote several functions:

• Rotation around X axis
• Rotation around Y axis
• Rotation around Z axis
• Rotation around all axes
• Rotation around any given axis
• Rotation from normal vector to normal vector

Apparently the 5th function is enough, because for example "Rotation around X axis" can be replace by rotation around (1,0,0), and "Rotation around all axes" is merely the product of 3 matrices. BUT, one should use the specific function he or she needs, because it is more efficient.
For matrix operations you can use this post.

Code:

public static Matrix4 GetRotationMatrixX(double angle)
{
    if (angle == 0.0)
    {
        return Matrix4.I;
    }
    float sin = (float)Math.Sin(angle);
    float cos = (float)Math.Cos(angle);
    return new Matrix4(new float[4, 4] {
        { 1.0f, 0.0f, 0.0f, 0.0f }, 
        { 0.0f, cos, -sin, 0.0f }, 
        { 0.0f, sin, cos, 0.0f }, 
        { 0.0f, 0.0f, 0.0f, 1.0f } });
}

public static Matrix4 GetRotationMatrixY(double angle)
{
    if (angle == 0.0)
    {
        return Matrix4.I;
    }
    float sin = (float)Math.Sin(angle);
    float cos = (float)Math.Cos(angle);
    return new Matrix4(new float[4, 4] {
        { cos, 0.0f, sin, 0.0f }, 
        { 0.0f, 1.0f, 0.0f, 0.0f }, 
        { -sin, 0.0f, cos, 0.0f }, 
        { 0.0f, 0.0f, 0.0f, 1.0f } });
}

public static Matrix4 GetRotationMatrixZ(double angle)
{
    if (angle == 0.0)
    {
        return Matrix4.I;
    }
    float sin = (float)Math.Sin(angle);
    float cos = (float)Math.Cos(angle);
    return new Matrix4(new float[4, 4] {
        { cos, -sin, 0.0f, 0.0f }, 
        { sin, cos, 0.0f, 0.0f }, 
        { 0.0f, 0.0f, 1.0f, 0.0f }, 
        { 0.0f, 0.0f, 0.0f, 1.0f } });
}

public static Matrix4 GetRotationMatrix(double ax, double ay, double az)
{
    Matrix4 my = null;
    Matrix4 mz = null;
    Matrix4 result = null;
    if (ax != 0.0)
    {
        result = GetRotationMatrixX(ax);
    }
    if (ay != 0.0)
    {
        my = GetRotationMatrixY(ay);
    }
    if (az != 0.0)
    {
        mz = GetRotationMatrixZ(az);
    }
    if (my != null)
    {
        if (result != null)
        {
            result *= my;
        }
        else
        {
            result = my;
        }
    }
    if (mz != null)
    {
        if (result != null)
        {
            result *= mz;
        }
        else
        {
            result = mz;
        }
    }
    if (result != null)
    {
        return result;
    }
    else
    {
        return Matrix4.I;
    }
}

public static Matrix4 GetRotationMatrix(Vector3 axis, double angle)
{
    if (angle == 0.0)
    {
        return Matrix4.I;
    }

    float x = axis.x;
    float y = axis.y;
    float z = axis.z;
    float sin = (float)Math.Sin(angle);
    float cos = (float)Math.Cos(angle);
    float xx = x * x;
    float yy = y * y;
    float zz = z * z;
    float xy = x * y;
    float xz = x * z;
    float yz = y * z;

    float[,] matrix = new float[4, 4];

    matrix[0, 0] = xx + (1 - xx) * cos;
    matrix[1, 0] = xy * (1 - cos) + z * sin;
    matrix[2, 0] = xz * (1 - cos) - y * sin;
    matrix[3, 0] = 0.0f;

    matrix[0, 1] = xy * (1 - cos) - z * sin;
    matrix[1, 1] = yy + (1 - yy) * cos;
    matrix[2, 1] = yz * (1 - cos) + x * sin;
    matrix[3, 1] = 0.0f;

    matrix[0, 2] = xz * (1 - cos) + y * sin;
    matrix[1, 2] = yz * (1 - cos) - x * sin;
    matrix[2, 2] = zz + (1 - zz) * cos;
    matrix[3, 2] = 0.0f;

    matrix[3, 0] = 0.0f;
    matrix[3, 1] = 0.0f;
    matrix[3, 2] = 0.0f;
    matrix[3, 3] = 1.0f;

    return new Matrix4(matrix);
}

/// <param name="source">Should be normalized</param>
/// <param name="destination">Should be normalized</param>
public static Matrix4 GetRotationMatrix(Vector3 source, Vector3 destination)
{
    Vector3 rotaxis = Vector3.CrossProduct(source, destination);
    if (rotaxis != Vector3.Zero)
    {
        rotaxis.Normalize();
        float cos = source.DotProduct(destination);
        double angle = Math.Acos(cos);
        return GetRotationMatrix(rotaxis, angle);
    }
    else
    {
        return Matrix4.I;
    }
}

Step-by-Step Math

Have you ever given up working on a math problem because you couldn’t figure out the next step? Wolfram|Alpha can guide you step by step through the process of solving many mathematical problems, from solving a simple quadratic equation to taking the integral of a complex function.

From here:
http://blog.wolframalpha.com/2009/12/01/step-by-step-math/

Example:
http://www.wolframalpha.com/input/?i=integrate+x+Sqrt[1-Sqrt[x]]