/src/libjxl/lib/jxl/base/matrix_ops.h
Line | Count | Source |
1 | | // Copyright (c) the JPEG XL Project Authors. All rights reserved. |
2 | | // |
3 | | // Use of this source code is governed by a BSD-style |
4 | | // license that can be found in the LICENSE file. |
5 | | |
6 | | #ifndef LIB_JXL_BASE_MATRIX_OPS_H_ |
7 | | #define LIB_JXL_BASE_MATRIX_OPS_H_ |
8 | | |
9 | | // 3x3 matrix operations. |
10 | | |
11 | | #include <array> |
12 | | #include <cmath> // abs |
13 | | #include <cstddef> |
14 | | |
15 | | #include "lib/jxl/base/status.h" |
16 | | |
17 | | namespace jxl { |
18 | | |
19 | | using Vector3 = std::array<float, 3>; |
20 | | using Vector3d = std::array<double, 3>; |
21 | | using Matrix3x3 = std::array<Vector3, 3>; |
22 | | using Matrix3x3d = std::array<Vector3d, 3>; |
23 | | |
24 | | // Computes C = A * B, where A, B, C are 3x3 matrices. |
25 | | template <typename Matrix> |
26 | 415k | void Mul3x3Matrix(const Matrix& a, const Matrix& b, Matrix& c) { |
27 | 1.66M | for (size_t x = 0; x < 3; x++) { |
28 | 1.24M | alignas(16) Vector3d temp{b[0][x], b[1][x], b[2][x]}; // transpose |
29 | 4.98M | for (size_t y = 0; y < 3; y++) { |
30 | 3.73M | c[y][x] = a[y][0] * temp[0] + a[y][1] * temp[1] + a[y][2] * temp[2]; |
31 | 3.73M | } |
32 | 1.24M | } |
33 | 415k | } |
34 | | |
35 | | // Computes C = A * B, where A is 3x3 matrix and B is vector. |
36 | | template <typename Matrix, typename Vector> |
37 | 343k | void Mul3x3Vector(const Matrix& a, const Vector& b, Vector& c) { |
38 | 1.37M | for (size_t y = 0; y < 3; y++) { |
39 | 1.02M | double e = 0; |
40 | 4.11M | for (size_t x = 0; x < 3; x++) { |
41 | 3.08M | e += static_cast<double>(a[y][x]) * b[x]; |
42 | 3.08M | } |
43 | 1.02M | c[y] = e; |
44 | 1.02M | } |
45 | 343k | } |
46 | | |
47 | | // Inverts a 3x3 matrix in place. |
48 | | template <typename Matrix> |
49 | 72.2k | Status Inv3x3Matrix(Matrix& matrix) { |
50 | | // Intermediate computation is done in double precision. |
51 | 72.2k | Matrix3x3d temp; |
52 | 72.2k | temp[0][0] = static_cast<double>(matrix[1][1]) * matrix[2][2] - |
53 | 72.2k | static_cast<double>(matrix[1][2]) * matrix[2][1]; |
54 | 72.2k | temp[0][1] = static_cast<double>(matrix[0][2]) * matrix[2][1] - |
55 | 72.2k | static_cast<double>(matrix[0][1]) * matrix[2][2]; |
56 | 72.2k | temp[0][2] = static_cast<double>(matrix[0][1]) * matrix[1][2] - |
57 | 72.2k | static_cast<double>(matrix[0][2]) * matrix[1][1]; |
58 | 72.2k | temp[1][0] = static_cast<double>(matrix[1][2]) * matrix[2][0] - |
59 | 72.2k | static_cast<double>(matrix[1][0]) * matrix[2][2]; |
60 | 72.2k | temp[1][1] = static_cast<double>(matrix[0][0]) * matrix[2][2] - |
61 | 72.2k | static_cast<double>(matrix[0][2]) * matrix[2][0]; |
62 | 72.2k | temp[1][2] = static_cast<double>(matrix[0][2]) * matrix[1][0] - |
63 | 72.2k | static_cast<double>(matrix[0][0]) * matrix[1][2]; |
64 | 72.2k | temp[2][0] = static_cast<double>(matrix[1][0]) * matrix[2][1] - |
65 | 72.2k | static_cast<double>(matrix[1][1]) * matrix[2][0]; |
66 | 72.2k | temp[2][1] = static_cast<double>(matrix[0][1]) * matrix[2][0] - |
67 | 72.2k | static_cast<double>(matrix[0][0]) * matrix[2][1]; |
68 | 72.2k | temp[2][2] = static_cast<double>(matrix[0][0]) * matrix[1][1] - |
69 | 72.2k | static_cast<double>(matrix[0][1]) * matrix[1][0]; |
70 | 72.2k | double det = matrix[0][0] * temp[0][0] + matrix[0][1] * temp[1][0] + |
71 | 72.2k | matrix[0][2] * temp[2][0]; |
72 | 72.2k | if (std::abs(det) < 1e-10) { |
73 | 0 | return JXL_FAILURE("Matrix determinant is too close to 0"); |
74 | 0 | } |
75 | 72.2k | double idet = 1.0 / det; |
76 | 288k | for (size_t j = 0; j < 3; j++) { |
77 | 866k | for (size_t i = 0; i < 3; i++) { |
78 | 650k | matrix[j][i] = temp[j][i] * idet; |
79 | 650k | } |
80 | 216k | } |
81 | 72.2k | return true; |
82 | 72.2k | } jxl::Status jxl::Inv3x3Matrix<std::__1::array<std::__1::array<float, 3ul>, 3ul> >(std::__1::array<std::__1::array<float, 3ul>, 3ul>&) Line | Count | Source | 49 | 71.9k | Status Inv3x3Matrix(Matrix& matrix) { | 50 | | // Intermediate computation is done in double precision. | 51 | 71.9k | Matrix3x3d temp; | 52 | 71.9k | temp[0][0] = static_cast<double>(matrix[1][1]) * matrix[2][2] - | 53 | 71.9k | static_cast<double>(matrix[1][2]) * matrix[2][1]; | 54 | 71.9k | temp[0][1] = static_cast<double>(matrix[0][2]) * matrix[2][1] - | 55 | 71.9k | static_cast<double>(matrix[0][1]) * matrix[2][2]; | 56 | 71.9k | temp[0][2] = static_cast<double>(matrix[0][1]) * matrix[1][2] - | 57 | 71.9k | static_cast<double>(matrix[0][2]) * matrix[1][1]; | 58 | 71.9k | temp[1][0] = static_cast<double>(matrix[1][2]) * matrix[2][0] - | 59 | 71.9k | static_cast<double>(matrix[1][0]) * matrix[2][2]; | 60 | 71.9k | temp[1][1] = static_cast<double>(matrix[0][0]) * matrix[2][2] - | 61 | 71.9k | static_cast<double>(matrix[0][2]) * matrix[2][0]; | 62 | 71.9k | temp[1][2] = static_cast<double>(matrix[0][2]) * matrix[1][0] - | 63 | 71.9k | static_cast<double>(matrix[0][0]) * matrix[1][2]; | 64 | 71.9k | temp[2][0] = static_cast<double>(matrix[1][0]) * matrix[2][1] - | 65 | 71.9k | static_cast<double>(matrix[1][1]) * matrix[2][0]; | 66 | 71.9k | temp[2][1] = static_cast<double>(matrix[0][1]) * matrix[2][0] - | 67 | 71.9k | static_cast<double>(matrix[0][0]) * matrix[2][1]; | 68 | 71.9k | temp[2][2] = static_cast<double>(matrix[0][0]) * matrix[1][1] - | 69 | 71.9k | static_cast<double>(matrix[0][1]) * matrix[1][0]; | 70 | 71.9k | double det = matrix[0][0] * temp[0][0] + matrix[0][1] * temp[1][0] + | 71 | 71.9k | matrix[0][2] * temp[2][0]; | 72 | 71.9k | if (std::abs(det) < 1e-10) { | 73 | 0 | return JXL_FAILURE("Matrix determinant is too close to 0"); | 74 | 0 | } | 75 | 71.9k | double idet = 1.0 / det; | 76 | 287k | for (size_t j = 0; j < 3; j++) { | 77 | 863k | for (size_t i = 0; i < 3; i++) { | 78 | 647k | matrix[j][i] = temp[j][i] * idet; | 79 | 647k | } | 80 | 215k | } | 81 | 71.9k | return true; | 82 | 71.9k | } |
jxl::Status jxl::Inv3x3Matrix<std::__1::array<std::__1::array<double, 3ul>, 3ul> >(std::__1::array<std::__1::array<double, 3ul>, 3ul>&) Line | Count | Source | 49 | 258 | Status Inv3x3Matrix(Matrix& matrix) { | 50 | | // Intermediate computation is done in double precision. | 51 | 258 | Matrix3x3d temp; | 52 | 258 | temp[0][0] = static_cast<double>(matrix[1][1]) * matrix[2][2] - | 53 | 258 | static_cast<double>(matrix[1][2]) * matrix[2][1]; | 54 | 258 | temp[0][1] = static_cast<double>(matrix[0][2]) * matrix[2][1] - | 55 | 258 | static_cast<double>(matrix[0][1]) * matrix[2][2]; | 56 | 258 | temp[0][2] = static_cast<double>(matrix[0][1]) * matrix[1][2] - | 57 | 258 | static_cast<double>(matrix[0][2]) * matrix[1][1]; | 58 | 258 | temp[1][0] = static_cast<double>(matrix[1][2]) * matrix[2][0] - | 59 | 258 | static_cast<double>(matrix[1][0]) * matrix[2][2]; | 60 | 258 | temp[1][1] = static_cast<double>(matrix[0][0]) * matrix[2][2] - | 61 | 258 | static_cast<double>(matrix[0][2]) * matrix[2][0]; | 62 | 258 | temp[1][2] = static_cast<double>(matrix[0][2]) * matrix[1][0] - | 63 | 258 | static_cast<double>(matrix[0][0]) * matrix[1][2]; | 64 | 258 | temp[2][0] = static_cast<double>(matrix[1][0]) * matrix[2][1] - | 65 | 258 | static_cast<double>(matrix[1][1]) * matrix[2][0]; | 66 | 258 | temp[2][1] = static_cast<double>(matrix[0][1]) * matrix[2][0] - | 67 | 258 | static_cast<double>(matrix[0][0]) * matrix[2][1]; | 68 | 258 | temp[2][2] = static_cast<double>(matrix[0][0]) * matrix[1][1] - | 69 | 258 | static_cast<double>(matrix[0][1]) * matrix[1][0]; | 70 | 258 | double det = matrix[0][0] * temp[0][0] + matrix[0][1] * temp[1][0] + | 71 | 258 | matrix[0][2] * temp[2][0]; | 72 | 258 | if (std::abs(det) < 1e-10) { | 73 | 0 | return JXL_FAILURE("Matrix determinant is too close to 0"); | 74 | 0 | } | 75 | 258 | double idet = 1.0 / det; | 76 | 1.03k | for (size_t j = 0; j < 3; j++) { | 77 | 3.09k | for (size_t i = 0; i < 3; i++) { | 78 | 2.32k | matrix[j][i] = temp[j][i] * idet; | 79 | 2.32k | } | 80 | 774 | } | 81 | 258 | return true; | 82 | 258 | } |
|
83 | | |
84 | | } // namespace jxl |
85 | | |
86 | | #endif // LIB_JXL_BASE_MATRIX_OPS_H_ |