enc_linalg.cc (1100B)
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 #include "lib/jxl/enc_linalg.h" 7 8 #include <cmath> 9 10 #include "lib/jxl/base/status.h" 11 12 namespace jxl { 13 14 void ConvertToDiagonal(const Matrix2x2& A, Vector2& diag, Matrix2x2& U) { 15 // Check A is symmetric. 16 JXL_DASSERT(std::abs(A[0][1] - A[1][0]) < 1e-15); 17 18 double b = -(A[0][0] + A[1][1]); 19 double c = A[0][0] * A[1][1] - A[0][1] * A[0][1]; 20 double d = b * b - 4.0 * c; 21 22 if (std::abs(A[0][1]) < 1e-10 || d < 0) { 23 // Already diagonal. 24 diag[0] = A[0][0]; 25 diag[1] = A[1][1]; 26 U[0][0] = U[1][1] = 1.0; 27 U[0][1] = U[1][0] = 0.0; 28 return; 29 } 30 31 double sqd = std::sqrt(d); 32 double l1 = (-b - sqd) * 0.5; 33 double l2 = (-b + sqd) * 0.5; 34 35 Vector2 v1 = {A[0][0] - l1, A[1][0]}; 36 double v1n = 1.0 / std::hypot(v1[0], v1[1]); 37 v1[0] = v1[0] * v1n; 38 v1[1] = v1[1] * v1n; 39 40 diag[0] = l1; 41 diag[1] = l2; 42 43 U[0][0] = v1[1]; 44 U[0][1] = -v1[0]; 45 U[1][0] = v1[0]; 46 U[1][1] = v1[1]; 47 } 48 49 } // namespace jxl