Coverage Report

Created: 2025-08-12 07:37

/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
70.6M
void Mul3x3Matrix(const Matrix& a, const Matrix& b, Matrix& c) {
27
282M
  for (size_t x = 0; x < 3; x++) {
28
211M
    alignas(16) Vector3d temp{b[0][x], b[1][x], b[2][x]};  // transpose
29
847M
    for (size_t y = 0; y < 3; y++) {
30
635M
      c[y][x] = a[y][0] * temp[0] + a[y][1] * temp[1] + a[y][2] * temp[2];
31
635M
    }
32
211M
  }
33
70.6M
}
34
35
// Computes C = A * B, where A is 3x3 matrix and B is vector.
36
template <typename Matrix, typename Vector>
37
53.0M
void Mul3x3Vector(const Matrix& a, const Vector& b, Vector& c) {
38
212M
  for (size_t y = 0; y < 3; y++) {
39
159M
    double e = 0;
40
637M
    for (size_t x = 0; x < 3; x++) {
41
477M
      e += a[y][x] * b[x];
42
477M
    }
43
159M
    c[y] = e;
44
159M
  }
45
53.0M
}
46
47
// Inverts a 3x3 matrix in place.
48
template <typename Matrix>
49
17.5M
Status Inv3x3Matrix(Matrix& matrix) {
50
  // Intermediate computation is done in double precision.
51
17.5M
  Matrix3x3d temp;
52
17.5M
  temp[0][0] = static_cast<double>(matrix[1][1]) * matrix[2][2] -
53
17.5M
               static_cast<double>(matrix[1][2]) * matrix[2][1];
54
17.5M
  temp[0][1] = static_cast<double>(matrix[0][2]) * matrix[2][1] -
55
17.5M
               static_cast<double>(matrix[0][1]) * matrix[2][2];
56
17.5M
  temp[0][2] = static_cast<double>(matrix[0][1]) * matrix[1][2] -
57
17.5M
               static_cast<double>(matrix[0][2]) * matrix[1][1];
58
17.5M
  temp[1][0] = static_cast<double>(matrix[1][2]) * matrix[2][0] -
59
17.5M
               static_cast<double>(matrix[1][0]) * matrix[2][2];
60
17.5M
  temp[1][1] = static_cast<double>(matrix[0][0]) * matrix[2][2] -
61
17.5M
               static_cast<double>(matrix[0][2]) * matrix[2][0];
62
17.5M
  temp[1][2] = static_cast<double>(matrix[0][2]) * matrix[1][0] -
63
17.5M
               static_cast<double>(matrix[0][0]) * matrix[1][2];
64
17.5M
  temp[2][0] = static_cast<double>(matrix[1][0]) * matrix[2][1] -
65
17.5M
               static_cast<double>(matrix[1][1]) * matrix[2][0];
66
17.5M
  temp[2][1] = static_cast<double>(matrix[0][1]) * matrix[2][0] -
67
17.5M
               static_cast<double>(matrix[0][0]) * matrix[2][1];
68
17.5M
  temp[2][2] = static_cast<double>(matrix[0][0]) * matrix[1][1] -
69
17.5M
               static_cast<double>(matrix[0][1]) * matrix[1][0];
70
17.5M
  double det = matrix[0][0] * temp[0][0] + matrix[0][1] * temp[1][0] +
71
17.5M
               matrix[0][2] * temp[2][0];
72
17.5M
  if (std::abs(det) < 1e-10) {
73
3
    return JXL_FAILURE("Matrix determinant is too close to 0");
74
3
  }
75
17.5M
  double idet = 1.0 / det;
76
70.3M
  for (size_t j = 0; j < 3; j++) {
77
210M
    for (size_t i = 0; i < 3; i++) {
78
158M
      matrix[j][i] = temp[j][i] * idet;
79
158M
    }
80
52.7M
  }
81
17.5M
  return true;
82
17.5M
}
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
17.5M
Status Inv3x3Matrix(Matrix& matrix) {
50
  // Intermediate computation is done in double precision.
51
17.5M
  Matrix3x3d temp;
52
17.5M
  temp[0][0] = static_cast<double>(matrix[1][1]) * matrix[2][2] -
53
17.5M
               static_cast<double>(matrix[1][2]) * matrix[2][1];
54
17.5M
  temp[0][1] = static_cast<double>(matrix[0][2]) * matrix[2][1] -
55
17.5M
               static_cast<double>(matrix[0][1]) * matrix[2][2];
56
17.5M
  temp[0][2] = static_cast<double>(matrix[0][1]) * matrix[1][2] -
57
17.5M
               static_cast<double>(matrix[0][2]) * matrix[1][1];
58
17.5M
  temp[1][0] = static_cast<double>(matrix[1][2]) * matrix[2][0] -
59
17.5M
               static_cast<double>(matrix[1][0]) * matrix[2][2];
60
17.5M
  temp[1][1] = static_cast<double>(matrix[0][0]) * matrix[2][2] -
61
17.5M
               static_cast<double>(matrix[0][2]) * matrix[2][0];
62
17.5M
  temp[1][2] = static_cast<double>(matrix[0][2]) * matrix[1][0] -
63
17.5M
               static_cast<double>(matrix[0][0]) * matrix[1][2];
64
17.5M
  temp[2][0] = static_cast<double>(matrix[1][0]) * matrix[2][1] -
65
17.5M
               static_cast<double>(matrix[1][1]) * matrix[2][0];
66
17.5M
  temp[2][1] = static_cast<double>(matrix[0][1]) * matrix[2][0] -
67
17.5M
               static_cast<double>(matrix[0][0]) * matrix[2][1];
68
17.5M
  temp[2][2] = static_cast<double>(matrix[0][0]) * matrix[1][1] -
69
17.5M
               static_cast<double>(matrix[0][1]) * matrix[1][0];
70
17.5M
  double det = matrix[0][0] * temp[0][0] + matrix[0][1] * temp[1][0] +
71
17.5M
               matrix[0][2] * temp[2][0];
72
17.5M
  if (std::abs(det) < 1e-10) {
73
3
    return JXL_FAILURE("Matrix determinant is too close to 0");
74
3
  }
75
17.5M
  double idet = 1.0 / det;
76
70.3M
  for (size_t j = 0; j < 3; j++) {
77
210M
    for (size_t i = 0; i < 3; i++) {
78
158M
      matrix[j][i] = temp[j][i] * idet;
79
158M
    }
80
52.7M
  }
81
17.5M
  return true;
82
17.5M
}
Unexecuted instantiation: jxl::Status jxl::Inv3x3Matrix<std::__1::array<std::__1::array<double, 3ul>, 3ul> >(std::__1::array<std::__1::array<double, 3ul>, 3ul>&)
83
84
}  // namespace jxl
85
86
#endif  // LIB_JXL_BASE_MATRIX_OPS_H_