vb01/matrix.cpp

#include "matrix.h"

#include <iostream>

using namespace std;

namespace vb01{
	Matrix::Matrix(float **matrix, int rows, int columns){
		this->dim[0] = rows;
		this->dim[1] = columns;

		this->matrix = new float*[rows];
		for(int i = 0; i < dim[0]; i++){
			this->matrix[i] = new float[columns];
			for(int j = 0; j < dim[1]; j++)
				this->matrix[i][j] = matrix[i][j];
		}
	}

	Matrix::Matrix(Vector3 v1, Vector3 v2, Vector3 v3){
		dim[0] = 3;
		dim[1] = 3;

		matrix=new float*[3];

		matrix[0] = new float[3];
		matrix[0][0] = v1.x;
		matrix[0][1] = v1.y;
		matrix[0][2] = v1.z;
		matrix[1] = new float[3];
		matrix[1][0] = v2.x;
		matrix[1][1] = v2.y;
		matrix[1][2] = v2.z;
		matrix[2] = new float[3];
		matrix[2][0] = v3.x;
		matrix[2][1] = v3.y;
		matrix[2][2] = v3.z;
	}

	Matrix& Matrix::operator+(Matrix &m){
		if(dim[0] == m.dim[0] && dim[1] == m.dim[1]){
			for(int i = 0; i < dim[0]; i++)
				for(int j = 0; j < dim[1]; j++)
					matrix[i][j] += m.matrix[i][j];
			return *this;
		}
		else{
			cout << "Incompatible matrix dimensions for addition.\n";
			exit(-1);
		}
	}

	Matrix& Matrix::operator-(Matrix &m){
		if(dim[0] == m.dim[0] && dim[1] == m.dim[1]){
			for(int i = 0; i < dim[0]; i++)
				for(int j = 0; j < dim[1]; j++)
					matrix[i][j] -= m.matrix[i][j];
			return *this;
		}
		else{
			cout << "Incompatible matrix dimensions for subtraction.\n";
			exit(-1);
		}
	}

	Matrix Matrix::operator*(Matrix &m){
		if(dim[1] == m.dim[0]){
			float **matrix = new float*[dim[0]];
			for(int i = 0; i < dim[0]; i++){
				matrix[i] = new float[m.dim[1]];
				for(int j = 0;j < m.dim[1]; j++)
					matrix[i][j] = 0.;
			}

			for(int i = 0; i < dim[0]; i++)
				for(int j = 0; j < dim[1]; j++)
					for(int k = 0; k < dim[0]; k++)
						matrix[i][j] += this->matrix[i][k] * m.matrix[k][j];
			return Matrix(matrix, dim[0], dim[1]);
		}
		else{
			cout << "Incompatible matrix dimensions for multiplication.\n";
			exit(-1);
		}
	}

	Vector3 Matrix::operator* (Vector3 v){
		float **matrix = new float*[3];
		matrix[0] = new float[1];
		matrix[0][0] = v.x;
		matrix[1] = new float[1];
		matrix[1][0] = v.y;
		matrix[2] = new float[1];
		matrix[2][0] = v.z;

		Matrix vMat(matrix, 3, 1);
		vMat = (*this) * vMat;
		return Vector3(vMat[0][0], vMat[1][0], vMat[2][0]);
	}

	template<typename T> Matrix& Matrix::operator+(T s){
		for(int i = 0; i < dim[0]; i++)
			for(int j = 0; j < dim[1]; j++)
				matrix[i][j] += s;
		return *this;
	}

	template<typename T> Matrix& Matrix::operator-(T s){
		for(int i = 0; i < dim[0]; i++)
			for(int j = 0; j < dim[1]; j++)
				matrix[i][j] -= s;
		return *this;
	}

	template<typename T> Matrix& Matrix::operator*(T s){
		for(int i = 0; i < dim[0]; i++)
			for(int j = 0; j < dim[1]; j++)
				matrix[i][j] *= s;
		return *this;
	}

	template<typename T> Matrix& Matrix::operator/(T s){
		for(int i = 0; i < dim[0]; i++)
			for(int j = 0; j < dim[1]; j++)
				matrix[i][j] /= s;
		return *this;
	}

	void Matrix::transpose(){
		float **matrix = new float*[dim[1]];
		for(int i = 0; i < dim[1]; i++)
			matrix[i] = new float[dim[0]];
		for(int i = 0; i < dim[0]; i++)
			for(int j = 0; j < dim[1]; j++)
				matrix[i][j] = this->matrix[j][i];
		this->matrix = matrix;
	}

	Matrix Matrix::getSubMatrix(int row, int column){
		int newRows = dim[0] - 1, newColumns = dim[1] - 1;
		float **matrix = new float*[dim[0]];
		bool offsetRows = false, offsetColumns = false;

		for(int i = 0; i < dim[1]; i++){
			matrix[i] = new float[dim[1]];

			for(int j = 0; j < dim[1]; j++)
				if(i != row && j != column){
					int a = (offsetRows && i != 0 ? i - 1 : i), b = (offsetColumns && j != 0 ? j - 1 : j);
					matrix[a][b] = this->matrix[i][j];
				}

				else if(i == row) offsetRows = true;
				else if(j == column) offsetColumns = true;
		}
		return Matrix(matrix, newRows, newColumns);
	}

	void Matrix::invert(){
		transpose();
		float **matrix = new float*[dim[0]];

		for(int i = 0; i < dim[0]; i++){
			matrix[i] = new float[dim[1]];
			for(int j = 0; j < dim[1]; j++)
				matrix[i][j] = getCofactor(i, j);
		}

		for(int i = 0; i < dim[0]; i++)
			delete this->matrix[i];
		delete[] this->matrix;

		this->matrix = matrix;
		*this / getDeterminant();
	}

	float Matrix::getMinor(int row, int column){
		return getSubMatrix(row, column).getDeterminant();
	}

	float Matrix::getCofactor(int row, int column){
		return ((row + column + 2) % 2 == 0 ? 1 : -1) * getMinor(row, column);	
	}

	float Matrix::getDeterminant(){
		if(dim[0] != dim[1]){
			cout << "Matrix not square.\n";
			exit(-1);
		}
		if(dim[0] == 1)
			return matrix[0][0];
		else if(dim[0] == 2)
			return matrix[0][0] * matrix[1][1] - matrix[1][0] * matrix[0][1];
		else{
			float det = 0;
			for(int i = 0; i < dim[0]; i++)
				det += matrix[i][0] * getCofactor(i, 0);
			return det;
		}
	}
}