godot/thirdparty/manifold/src/svd.h

// MIT License

// Copyright (c) 2019 wi-re
// Copyright 2023 The Manifold Authors.

// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:

// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.

// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.

// Modified from https://github.com/wi-re/tbtSVD, removing CUDA dependence and
// approximate inverse square roots.

#include <cmath>

#include "manifold/common.h"

namespace {
mat3;
vec3;
vec4;

// Constants used for calculation of Givens quaternions
inline constexpr double _gamma = …;    // sqrt(8)+3;
inline constexpr double _cStar = …;  // cos(pi/8)
inline constexpr double _sStar = …;  // sin(pi/8)
// Threshold value
inline constexpr double _SVD_EPSILON = …;
// Iteration counts for Jacobi Eigen Analysis, influences precision
inline constexpr int JACOBI_STEPS = …;

// Helper function used to swap X with Y and Y with  X if c == true
inline void CondSwap(bool c, double& X, double& Y) { … }
// Helper function used to swap X with Y and Y with -X if c == true
inline void CondNegSwap(bool c, double& X, double& Y) { … }
// A simple symmetric 3x3 Matrix class (contains no storage for (0, 1) (0, 2)
// and (1, 2)
struct Symmetric3x3 { … };
// Helper struct to store 2 doubles to avoid OUT parameters on functions
struct Givens { … };
// Helper struct to store 2 Matrices to avoid OUT parameters on functions
struct QR { … };
// Calculates the squared norm of the vector.
inline double Dist2(vec3 v) { … }
// For an explanation of the math see
// http://pages.cs.wisc.edu/~sifakis/papers/SVD_TR1690.pdf Computing the
// Singular Value Decomposition of 3 x 3 matrices with minimal branching and
// elementary floating point operations See Algorithm 2 in reference. Given a
// matrix A this function returns the Givens quaternion (x and w component, y
// and z are 0)
inline Givens ApproximateGivensQuaternion(Symmetric3x3& A) { … }
// Function used to apply a Givens rotation S. Calculates the weights and
// updates the quaternion to contain the cumulative rotation
inline void JacobiConjugation(const int32_t x, const int32_t y, const int32_t z,
                              Symmetric3x3& S, vec4& q) { … }
// Function used to contain the Givens permutations and the loop of the jacobi
// steps controlled by JACOBI_STEPS Returns the quaternion q containing the
// cumulative result used to reconstruct S
inline mat3 JacobiEigenAnalysis(Symmetric3x3 S) { … }
// Implementation of Algorithm 3
inline void SortSingularValues(mat3& B, mat3& V) { … }
// Implementation of Algorithm 4
inline Givens QRGivensQuaternion(double a1, double a2) { … }
// Implements a QR decomposition of a Matrix, see Sec 4.2
inline QR QRDecomposition(mat3& B) { … }
}  // namespace

namespace manifold {

/**
 * The three matrices of a Singular Value Decomposition.
 */
struct SVDSet { … };

/**
 * Returns the Singular Value Decomposition of A: A = U * S * la::transpose(V).
 *
 * @param A The matrix to decompose.
 */
inline SVDSet SVD(mat3 A) { … }

/**
 * Returns the largest singular value of A.
 *
 * @param A The matrix to measure.
 */
inline double SpectralNorm(mat3 A) { … }
}  // namespace manifold