1
0
forked from josh/j3ml
Files
j3ml-fork/src/J3ML/LinearAlgebra/Matrix4x4.cpp
2024-02-01 17:23:48 -05:00

422 lines
18 KiB
C++

#include <J3ML/LinearAlgebra/Matrix4x4.h>
#include <J3ML/LinearAlgebra/Vector4.h>
namespace LinearAlgebra {
const Matrix4x4 Matrix4x4::Zero = Matrix4x4(0);
const Matrix4x4 Matrix4x4::Identity = Matrix4x4({1,0,0,0}, {0,1,0,0}, {0,0,1,0}, {0,0,0,1});
const Matrix4x4 Matrix4x4::NaN = Matrix4x4(NAN);
Matrix4x4::Matrix4x4(const Vector4 &r1, const Vector4 &r2, const Vector4 &r3, const Vector4 &r4) {
this->elems[0][0] = r1.x;
this->elems[0][1] = r1.y;
this->elems[0][2] = r1.z;
this->elems[0][3] = r1.w;
this->elems[1][0] = r2.x;
this->elems[1][1] = r2.y;
this->elems[1][2] = r2.z;
this->elems[1][3] = r2.w;
this->elems[2][0] = r3.x;
this->elems[2][1] = r3.y;
this->elems[2][2] = r3.z;
this->elems[2][3] = r3.w;
this->elems[3][0] = r4.x;
this->elems[3][1] = r4.y;
this->elems[3][2] = r4.z;
this->elems[3][3] = r4.w;
}
Matrix4x4::Matrix4x4(float val) {
this->elems[0][0] = val;
this->elems[0][1] = val;
this->elems[0][2] = val;
this->elems[0][3] = val;
this->elems[1][0] = val;
this->elems[1][1] = val;
this->elems[1][2] = val;
this->elems[1][3] = val;
this->elems[2][0] = val;
this->elems[2][1] = val;
this->elems[2][2] = val;
this->elems[2][3] = val;
this->elems[3][0] = val;
this->elems[3][1] = val;
this->elems[3][2] = val;
this->elems[3][3] = val;
}
Matrix4x4::Matrix4x4(float m00, float m01, float m02, float m03, float m10, float m11, float m12, float m13,
float m20, float m21, float m22, float m23, float m30, float m31, float m32, float m33) {
this->elems[0][0] = m00;
this->elems[0][1] = m01;
this->elems[0][2] = m02;
this->elems[0][3] = m03;
this->elems[1][0] = m10;
this->elems[1][1] = m11;
this->elems[1][2] = m12;
this->elems[1][3] = m13;
this->elems[2][0] = m20;
this->elems[2][1] = m21;
this->elems[2][2] = m22;
this->elems[2][3] = m23;
this->elems[3][0] = m30;
this->elems[3][1] = m31;
this->elems[3][2] = m32;
this->elems[3][3] = m33;
}
Matrix4x4::Matrix4x4(const Quaternion &orientation) {
}
void Matrix4x4::SetTranslatePart(float translateX, float translateY, float translateZ) {
elems[0][3] = translateX;
elems[1][3] = translateY;
elems[2][3] = translateZ;
}
void Matrix4x4::SetTranslatePart(const Vector3 &offset) {
elems[0][3] = offset.x;
elems[1][3] = offset.y;
elems[2][3] = offset.z;
}
void Matrix4x4::SetRotatePart(const Quaternion& q)
{
// See e.g. http://www.geometrictools.com/Documentation/LinearAlgebraicQuaternions.pdf .
const float x = q.x;
const float y = q.y;
const float z = q.z;
const float w = q.w;
elems[0][0] = 1 - 2*(y*y + z*z); elems[0][1] = 2*(x*y - z*w); elems[0][2] = 2*(x*z + y*w);
elems[1][0] = 2*(x*y + z*w); elems[1][1] = 1 - 2*(x*x + z*z); elems[1][2] = 2*(y*z - x*w);
elems[2][0] = 2*(x*z - y*w); elems[2][1] = 2*(y*z + x*w); elems[2][2] = 1 - 2*(x*x + y*y);
}
void Matrix4x4::Set3x3Part(const Matrix3x3& r)
{
At(0, 0) = r[0][0]; At(0, 1) = r[0][1]; At(0, 2) = r[0][2];
At(1, 0) = r[1][0]; At(1, 1) = r[1][1]; At(1, 2) = r[1][2];
At(2, 0) = r[2][0]; At(2, 1) = r[2][1]; At(2, 2) = r[2][2];
}
void Matrix4x4::SetRow(int row, const Vector3 &rowVector, float m_r3) {
SetRow(row, rowVector.x, rowVector.y, rowVector.z, m_r3);
}
void Matrix4x4::SetRow(int row, const Vector4 &rowVector) {
SetRow(row, rowVector.x, rowVector.y, rowVector.z, rowVector.w);
}
void Matrix4x4::SetRow(int row, float m_r0, float m_r1, float m_r2, float m_r3) {
elems[row][0] = m_r0;
elems[row][1] = m_r1;
elems[row][2] = m_r2;
elems[row][3] = m_r3;
}
Matrix4x4::Matrix4x4(const Quaternion &orientation, const Vector3 &translation) {
SetRotatePart(orientation);
SetTranslatePart(translation);
SetRow(3, 0, 0, 0, 1);
}
Matrix4x4 Matrix4x4::D3DOrthoProjLH(float n, float f, float h, float v) {
float p00 = 2.f / h; float p01 = 0; float p02 = 0; float p03 = 0.f;
float p10 = 0; float p11 = 2.f / v; float p12 = 0; float p13 = 0.f;
float p20 = 0; float p21 = 0; float p22 = 1.f / (f-n); float p23 = n / (n-f);
float p30 = 0; float p31 = 0; float p32 = 0.f; float p33 = 1.f;
return {p00, p01, p02, p03, p10, p11, p12, p13, p20, p21, p22, p23, p30, p31, p32, p33};
}
/** This function generates an orthographic projection matrix that maps from
the Direct3D view space to the Direct3D normalized viewport space as follows:
In Direct3D view space, we assume that the camera is positioned at the origin (0,0,0).
The camera looks directly towards the positive Z axis (0,0,1).
The -X axis spans to the left of the screen, +X goes to the right.
-Y goes to the bottom of the screen, +Y goes to the top.
After the transformation, we're in the Direct3D normalized viewport space as follows:
(-1,-1,0) is the bottom-left corner of the viewport at the near plane.
(1,1,0) is the top-right corner of the viewport at the near plane.
(0,0,0) is the center point at the near plane.
Coordinates with z=1 are at the far plane.
Examples:
(0,0,n) maps to (0,0,0).
(0,0,f) maps to (0,0,1).
(-h/2, -v/2, n) maps to (-1, -1, 0).
(h/2, v/2, f) maps to (1, 1, 1).
*/
Matrix4x4 Matrix4x4::D3DOrthoProjRH(float n, float f, float h, float v)
{
// D3DOrthoProjLH and D3DOrthoProjRH differ from each other in that the third column is negated.
// This corresponds to LH = RH * In, where In is a diagonal matrix with elements [1 1 -1 1].
float p00 = 2.f / h; float p01 = 0; float p02 = 0; float p03 = 0.f;
float p10 = 0; float p11 = 2.f / v; float p12 = 0; float p13 = 0.f;
float p20 = 0; float p21 = 0; float p22 = 1.f / (n-f); float p23 = n / (n-f);
float p30 = 0; float p31 = 0; float p32 = 0.f; float p33 = 1.f;
}
float Matrix4x4::At(int x, int y) const {
return elems[x][y];
}
Matrix4x4 Matrix4x4::operator*(const Matrix4x4 &rhs) const {
float r00 = At(0, 0) * rhs.At(0, 0) + At(0, 1) * rhs.At(1, 0) + At(0, 2) * rhs.At(2, 0) + At(0, 3) * rhs.At(3, 0);
float r01 = At(0, 0) * rhs.At(0, 1) + At(0, 1) * rhs.At(1, 1) + At(0, 2) * rhs.At(2, 1) + At(0, 3) * rhs.At(3, 1);
float r02 = At(0, 0) * rhs.At(0, 2) + At(0, 1) * rhs.At(1, 2) + At(0, 2) * rhs.At(2, 2) + At(0, 3) * rhs.At(3, 2);
float r03 = At(0, 0) * rhs.At(0, 3) + At(0, 1) * rhs.At(1, 3) + At(0, 2) * rhs.At(2, 3) + At(0, 3) * rhs.At(3, 3);
float r10 = At(1, 0) * rhs.At(0, 0) + At(1, 1) * rhs.At(1, 0) + At(1, 2) * rhs.At(2, 0) + At(1, 3) * rhs.At(3, 0);
float r11 = At(1, 0) * rhs.At(0, 1) + At(1, 1) * rhs.At(1, 1) + At(1, 2) * rhs.At(2, 1) + At(1, 3) * rhs.At(3, 1);
float r12 = At(1, 0) * rhs.At(0, 2) + At(1, 1) * rhs.At(1, 2) + At(1, 2) * rhs.At(2, 2) + At(1, 3) * rhs.At(3, 2);
float r13 = At(1, 0) * rhs.At(0, 3) + At(1, 1) * rhs.At(1, 3) + At(1, 2) * rhs.At(2, 3) + At(1, 3) * rhs.At(3, 3);
float r20 = At(2, 0) * rhs.At(0, 0) + At(2, 1) * rhs.At(1, 0) + At(2, 2) * rhs.At(2, 0) + At(2, 3) * rhs.At(3, 0);
float r21 = At(2, 0) * rhs.At(0, 1) + At(2, 1) * rhs.At(1, 1) + At(2, 2) * rhs.At(2, 1) + At(2, 3) * rhs.At(3, 1);
float r22 = At(2, 0) * rhs.At(0, 2) + At(2, 1) * rhs.At(1, 2) + At(2, 2) * rhs.At(2, 2) + At(2, 3) * rhs.At(3, 2);
float r23 = At(2, 0) * rhs.At(0, 3) + At(2, 1) * rhs.At(1, 3) + At(2, 2) * rhs.At(2, 3) + At(2, 3) * rhs.At(3, 3);
float r30 = At(3, 0) * rhs.At(0, 0) + At(3, 1) * rhs.At(1, 0) + At(3, 2) * rhs.At(2, 0) + At(3, 3) * rhs.At(3, 0);
float r31 = At(3, 0) * rhs.At(0, 1) + At(3, 1) * rhs.At(1, 1) + At(3, 2) * rhs.At(2, 1) + At(3, 3) * rhs.At(3, 1);
float r32 = At(3, 0) * rhs.At(0, 2) + At(3, 1) * rhs.At(1, 2) + At(3, 2) * rhs.At(2, 2) + At(3, 3) * rhs.At(3, 2);
float r33 = At(3, 0) * rhs.At(0, 3) + At(3, 1) * rhs.At(1, 3) + At(3, 2) * rhs.At(2, 3) + At(3, 3) * rhs.At(3, 3);
return {r00,r01,r02,r03, r10, r11, r12, r13, r20,r21,r22,r23, r30,r31,r32,r33};
}
Vector4 Matrix4x4::operator[](int row) {
return Vector4{elems[row][0], elems[row][1], elems[row][2], elems[row][3]};
}
Matrix4x4 Matrix4x4::operator*(const Matrix3x3 &rhs) const {
float r00 = At(0, 0) * rhs.At(0, 0) + At(0, 1) * rhs.At(1, 0) + At(0, 2) * rhs.At(2, 0);
float r01 = At(0, 0) * rhs.At(0, 1) + At(0, 1) * rhs.At(1, 1) + At(0, 2) * rhs.At(2, 1);
float r02 = At(0, 0) * rhs.At(0, 2) + At(0, 1) * rhs.At(1, 2) + At(0, 2) * rhs.At(2, 2);
float r03 = At(0, 3);
float r10 = At(1, 0) * rhs.At(0, 0) + At(1, 1) * rhs.At(1, 0) + At(1, 2) * rhs.At(2, 0);
float r11 = At(1, 0) * rhs.At(0, 1) + At(1, 1) * rhs.At(1, 1) + At(1, 2) * rhs.At(2, 1);
float r12 = At(1, 0) * rhs.At(0, 2) + At(1, 1) * rhs.At(1, 2) + At(1, 2) * rhs.At(2, 2);
float r13 = At(1, 3);
float r20 = At(2, 0) * rhs.At(0, 0) + At(2, 1) * rhs.At(1, 0) + At(2, 2) * rhs.At(2, 0);
float r21 = At(2, 0) * rhs.At(0, 1) + At(2, 1) * rhs.At(1, 1) + At(2, 2) * rhs.At(2, 1);
float r22 = At(2, 0) * rhs.At(0, 2) + At(2, 1) * rhs.At(1, 2) + At(2, 2) * rhs.At(2, 2);
float r23 = At(2, 3);
float r30 = At(3, 0) * rhs.At(0, 0) + At(3, 1) * rhs.At(1, 0) + At(3, 2) * rhs.At(2, 0);
float r31 = At(3, 0) * rhs.At(0, 1) + At(3, 1) * rhs.At(1, 1) + At(3, 2) * rhs.At(2, 1);
float r32 = At(3, 0) * rhs.At(0, 2) + At(3, 1) * rhs.At(1, 2) + At(3, 2) * rhs.At(2, 2);
float r33 = At(3, 3);
return {r00,r01,r02,r03, r10, r11, r12, r13, r20,r21,r22,r23, r30,r31,r32,r33};
}
Matrix4x4 Matrix4x4::operator+() const { return *this; }
Matrix4x4 Matrix4x4::FromTranslation(const Vector3 &rhs) {
return Matrix4x4(1.f, 0, 0, rhs.x,
0, 1.f, 0, rhs.y,
0, 0, 1.f, rhs.z,
0, 0, 0, 1.f);
}
Matrix4x4 Matrix4x4::Translate(const Vector3 &rhs) const {
return *this * FromTranslation(rhs);
}
Vector3 Matrix4x4::Transform(const Vector3 &rhs) const {
return Transform(rhs.x, rhs.y, rhs.z);
}
Vector3 Matrix4x4::Transform(float tx, float ty, float tz) const {
return Vector3(At(0, 0) * tx + At(0, 1) * ty + At(0, 2) * tz + At(0, 3),
At(1, 0) * tx + At(1, 1) * ty + At(1, 2) * tz + At(1, 3),
At(2, 0) * tx + At(2, 1) * ty + At(2, 2) * tz + At(2, 3));
}
Vector2 Matrix4x4::Transform(float tx, float ty) const {
return Vector2(At(0, 0) * tx + At(0, 1) * ty + At(0, 2) + At(0, 3),
At(1, 0) * tx + At(1, 1) * ty + At(1, 2) + At(1, 3));
}
Vector2 Matrix4x4::Transform(const Vector2 &rhs) const {
return Transform(rhs.x, rhs.y);
}
Matrix4x4 &Matrix4x4::operator=(const Matrix3x3 &rhs) {
Set3x3Part(rhs);
SetTranslatePart(0,0,0);
SetRow(3, 0, 0, 0, 1);
return *this;
}
Matrix4x4 &Matrix4x4::operator=(const Quaternion &rhs) {
*this = rhs.ToMatrix4x4();
return *this;
}
float &Matrix4x4::At(int row, int col) {
return elems[row][col];
}
Matrix4x4 Matrix4x4::Inverse() const {
// Compute the inverse directly using Cramer's rule
// Warning: This method is numerically very unstable!
float d = Determinant();
d = 1.f / d;
float a11 = At(0, 0);float a12 = At(0, 1);float a13 = At(0, 2);float a14 = At(0, 3);
float a21 = At(1, 0);float a22 = At(1, 1);float a23 = At(1, 2);float a24 = At(1, 3);
float a31 = At(2, 0);float a32 = At(2, 1);float a33 = At(2, 2);float a34 = At(2, 3);
float a41 = At(3, 0);float a42 = At(3, 1);float a43 = At(3, 2);float a44 = At(3, 3);
Matrix4x4 i = {
d * (a22*a33*a44 + a23*a34*a42 + a24*a32*a43 - a22*a34*a43 - a23*a32*a44 - a24*a33*a42),
d * (a12*a34*a43 + a13*a32*a44 + a14*a33*a42 - a12*a33*a44 - a13*a34*a42 - a14*a32*a43),
d * (a12*a23*a44 + a13*a24*a42 + a14*a22*a43 - a12*a24*a43 - a13*a22*a44 - a14*a23*a42),
d * (a12*a24*a33 + a13*a22*a34 + a14*a23*a32 - a12*a23*a34 - a13*a24*a32 - a14*a22*a33),
d * (a21*a34*a43 + a23*a31*a44 + a24*a33*a41 - a21*a33*a44 - a23*a34*a41 - a24*a31*a43),
d * (a11*a33*a44 + a13*a34*a41 + a14*a31*a43 - a11*a34*a43 - a13*a31*a44 - a14*a33*a41),
d * (a11*a24*a43 + a13*a21*a44 + a14*a23*a41 - a11*a23*a44 - a13*a24*a41 - a14*a21*a43),
d * (a11*a23*a34 + a13*a24*a31 + a14*a21*a33 - a11*a24*a33 - a13*a21*a34 - a14*a23*a31),
d * (a21*a32*a44 + a22*a34*a41 + a24*a31*a42 - a21*a34*a42 - a22*a31*a44 - a24*a32*a41),
d * (a11*a34*a42 + a12*a31*a44 + a14*a32*a41 - a11*a32*a44 - a12*a34*a41 - a14*a31*a42),
d * (a11*a22*a44 + a12*a24*a41 + a14*a21*a42 - a11*a24*a42 - a12*a21*a44 - a14*a22*a41),
d * (a11*a24*a32 + a12*a21*a34 + a14*a22*a31 - a11*a22*a34 - a12*a24*a31 - a14*a21*a32),
d * (a21*a33*a42 + a22*a31*a43 + a23*a32*a41 - a21*a32*a43 - a22*a33*a41 - a23*a31*a42),
d * (a11*a32*a43 + a12*a33*a41 + a13*a31*a42 - a11*a33*a42 - a12*a31*a43 - a13*a32*a41),
d * (a11*a23*a42 + a12*a21*a43 + a13*a22*a41 - a11*a22*a43 - a12*a23*a41 - a13*a21*a42),
d * (a11*a22*a33 + a12*a23*a31 + a13*a21*a32 - a11*a23*a32 - a12*a21*a33 - a13*a22*a31)
};
return i;
}
float Matrix4x4::Minor(int i, int j) const {
int r0 = SKIPNUM(0, i);
int r1 = SKIPNUM(1, i);
int r2 = SKIPNUM(2, i);
int c0 = SKIPNUM(0, j);
int c1 = SKIPNUM(1, j);
int c2 = SKIPNUM(2, j);
float a = At(r0, c0);
float b = At(r0, c1);
float c = At(r0, c2);
float d = At(r1, c0);
float e = At(r1, c1);
float f = At(r1, c2);
float g = At(r2, c0);
float h = At(r2, c1);
float k = At(r2, c2);
return a*e*k + b*f*g + c*d*h - a*f*h - b*d*k - c*e*g;
}
float Matrix4x4::Determinant() const {
return At(0, 0) * Minor(0,0) - At(0, 1) * Minor(0,1) + At(0, 2) * Minor(0,2) - At(0, 3) * Minor(0,3);
}
float Matrix4x4::Determinant3x3() const {
const float a = elems[0][0];
const float b = elems[0][1];
const float c = elems[0][2];
const float d = elems[1][0];
const float e = elems[1][1];
const float f = elems[1][2];
const float g = elems[2][0];
const float h = elems[2][1];
const float i = elems[2][2];
return a*e*i + b*f*g + c*d*h - a*f*h - b*d*i - c*e*g;
}
Matrix3x3 Matrix4x4::GetRotatePart() const {
return Matrix3x3 {
At(0, 0), At(0, 1), At(0, 2),
At(1, 0), At(1, 1), At(1, 2),
At(2, 0), At(2, 1), At(2, 2)
};
}
Matrix4x4 Matrix4x4::Transpose() const {
Matrix4x4 copy;
copy.elems[0][0] = elems[0][0]; copy.elems[0][1] = elems[1][0]; copy.elems[0][2] = elems[2][0]; copy.elems[0][3] = elems[3][0];
copy.elems[1][0] = elems[0][1]; copy.elems[1][1] = elems[1][1]; copy.elems[1][2] = elems[2][1]; copy.elems[1][3] = elems[3][1];
copy.elems[2][0] = elems[0][2]; copy.elems[2][1] = elems[1][2]; copy.elems[2][2] = elems[2][2]; copy.elems[2][3] = elems[3][2];
copy.elems[3][0] = elems[0][3]; copy.elems[3][1] = elems[1][3]; copy.elems[3][2] = elems[2][3]; copy.elems[3][3] = elems[3][3];
return copy;
}
Vector4 Matrix4x4::Diagonal() const {
return Vector4{At(0, 0), At(1,1), At(2,2), At(3,3)};
}
Vector3 Matrix4x4::Diagonal3() const {
return Vector3 { At(0, 0), At(1,1), At(2,2) };
}
Vector3 Matrix4x4::WorldX() const {
return GetColumn3(0);
}
Vector3 Matrix4x4::WorldY() const {
return GetColumn3(1);
}
Vector3 Matrix4x4::WorldZ() const {
return GetColumn3(2);
}
bool Matrix4x4::IsFinite() const {
for(int iy = 0; iy < Rows; ++iy)
for(int ix = 0; ix < Cols; ++ix)
if (!std::isfinite(elems[iy][ix]))
return false;
return true;
}
Vector3 Matrix4x4::GetColumn3(int index) const {
return Vector3{At(0, index), At(1, index), At(2, index)};
}
Vector2 Matrix4x4::operator*(const Vector2 &rhs) const { return this->Transform(rhs);}
Vector3 Matrix4x4::operator*(const Vector3 &rhs) const { return this->Transform(rhs);}
Vector4 Matrix4x4::operator*(const Vector4 &rhs) const { return this->Transform(rhs);}
Vector4 Matrix4x4::Transform(float tx, float ty, float tz, float tw) const {
return Transform({tx, ty, tz, tw});
}
Vector4 Matrix4x4::Transform(const Vector4 &rhs) const {
return Vector4(At(0, 0) * rhs.x + At(0, 1) * rhs.y + At(0, 2) * rhs.z + At(0, 3) * rhs.w,
At(1, 0) * rhs.x + At(1, 1) * rhs.y + At(1, 2) * rhs.z + At(1, 3) * rhs.w,
At(2, 0) * rhs.x + At(2, 1) * rhs.y + At(2, 2) * rhs.z + At(2, 3) * rhs.w,
At(3, 0) * rhs.x + At(3, 1) * rhs.y + At(3, 2) * rhs.z + At(3, 3) * rhs.w);
}
Vector3 Matrix4x4::GetTranslatePart() const {
return GetColumn3(3);
}
}