Flipper/Applications/Official/source-OLDER/xMasterX/airmouse/tracking/util/matrixutils.cc

149 lines
4.6 KiB
C++

/*
* Copyright 2019 Google Inc. All Rights Reserved.
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#include "matrixutils.h"
#include "vectorutils.h"
namespace cardboard {
namespace {
// Returns true if the cofactor for a given row and column should be negated.
static bool IsCofactorNegated(int row, int col)
{
// Negated iff (row + col) is odd.
return ((row + col) & 1) != 0;
}
static double CofactorElement3(const Matrix3x3& m, int row, int col)
{
static const int index[3][2] = { { 1, 2 }, { 0, 2 }, { 0, 1 } };
const int i0 = index[row][0];
const int i1 = index[row][1];
const int j0 = index[col][0];
const int j1 = index[col][1];
const double cofactor = m(i0, j0) * m(i1, j1) - m(i0, j1) * m(i1, j0);
return IsCofactorNegated(row, col) ? -cofactor : cofactor;
}
// Multiplies a matrix and some type of column vector to
// produce another column vector of the same type.
Vector3 MultiplyMatrixAndVector(const Matrix3x3& m, const Vector3& v)
{
Vector3 result = Vector3::Zero();
for (int row = 0; row < 3; ++row) {
for (int col = 0; col < 3; ++col)
result[row] += m(row, col) * v[col];
}
return result;
}
// Sets the upper 3x3 of a Matrix to represent a 3D rotation.
void RotationMatrix3x3(const Rotation& r, Matrix3x3* matrix)
{
//
// Given a quaternion (a,b,c,d) where d is the scalar part, the 3x3 rotation
// matrix is:
//
// a^2 - b^2 - c^2 + d^2 2ab - 2cd 2ac + 2bd
// 2ab + 2cd -a^2 + b^2 - c^2 + d^2 2bc - 2ad
// 2ac - 2bd 2bc + 2ad -a^2 - b^2 + c^2 + d^2
//
const Vector<4>& quat = r.GetQuaternion();
const double aa = quat[0] * quat[0];
const double bb = quat[1] * quat[1];
const double cc = quat[2] * quat[2];
const double dd = quat[3] * quat[3];
const double ab = quat[0] * quat[1];
const double ac = quat[0] * quat[2];
const double bc = quat[1] * quat[2];
const double ad = quat[0] * quat[3];
const double bd = quat[1] * quat[3];
const double cd = quat[2] * quat[3];
Matrix3x3& m = *matrix;
m[0][0] = aa - bb - cc + dd;
m[0][1] = 2 * ab - 2 * cd;
m[0][2] = 2 * ac + 2 * bd;
m[1][0] = 2 * ab + 2 * cd;
m[1][1] = -aa + bb - cc + dd;
m[1][2] = 2 * bc - 2 * ad;
m[2][0] = 2 * ac - 2 * bd;
m[2][1] = 2 * bc + 2 * ad;
m[2][2] = -aa - bb + cc + dd;
}
} // anonymous namespace
Vector3 operator*(const Matrix3x3& m, const Vector3& v) { return MultiplyMatrixAndVector(m, v); }
Matrix3x3 CofactorMatrix(const Matrix3x3& m)
{
Matrix3x3 result;
for (int row = 0; row < 3; ++row) {
for (int col = 0; col < 3; ++col)
result(row, col) = CofactorElement3(m, row, col);
}
return result;
}
Matrix3x3 AdjugateWithDeterminant(const Matrix3x3& m, double* determinant)
{
const Matrix3x3 cofactor_matrix = CofactorMatrix(m);
if (determinant) {
*determinant = m(0, 0) * cofactor_matrix(0, 0) + m(0, 1) * cofactor_matrix(0, 1)
+ m(0, 2) * cofactor_matrix(0, 2);
}
return Transpose(cofactor_matrix);
}
// Returns the transpose of a matrix.
Matrix3x3 Transpose(const Matrix3x3& m)
{
Matrix3x3 result;
for (int row = 0; row < 3; ++row) {
for (int col = 0; col < 3; ++col)
result(row, col) = m(col, row);
}
return result;
}
Matrix3x3 InverseWithDeterminant(const Matrix3x3& m, double* determinant)
{
// The inverse is the adjugate divided by the determinant.
double det;
Matrix3x3 adjugate = AdjugateWithDeterminant(m, &det);
if (determinant)
*determinant = det;
if (det == 0)
return Matrix3x3::Zero();
else
return adjugate * (1 / det);
}
Matrix3x3 Inverse(const Matrix3x3& m) { return InverseWithDeterminant(m, nullptr); }
Matrix3x3 RotationMatrixNH(const Rotation& r)
{
Matrix3x3 m;
RotationMatrix3x3(r, &m);
return m;
}
} // namespace cardboard