# HG changeset patch # User Sebastien Jodogne # Date 1519990527 -3600 # Node ID 7105e51e49073c360c360c72893c5e3cf3ff2238 # Parent 62964bfbbc0024ad97b60f03e0a6fbf5db5d2c09 InvertMatrixUnsafe diff -r 62964bfbbc00 -r 7105e51e4907 Framework/Toolbox/LinearAlgebra.cpp --- a/Framework/Toolbox/LinearAlgebra.cpp Fri Mar 02 11:21:12 2018 +0100 +++ b/Framework/Toolbox/LinearAlgebra.cpp Fri Mar 02 12:35:27 2018 +0100 @@ -497,9 +497,9 @@ } } - - void InvertMatrix(Matrix& target, - const Matrix& source) + + bool InvertMatrixUnsafe(Matrix& target, + const Matrix& source) { if (source.size1() != source.size2()) { @@ -516,15 +516,14 @@ { // By convention, the inverse of the empty matrix, is itself the empty matrix target.resize(0, 0); - return; + return true; } double determinant = ComputeDeterminant(source); if (IsCloseToZero(determinant)) { - LOG(ERROR) << "Cannot invert singular matrix"; - throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); + return false; } double denominator = 1.0 / determinant; @@ -534,6 +533,8 @@ if (source.size1() == 1) { target(0, 0) = denominator; + + return true; } else if (source.size1() == 2) { @@ -542,6 +543,8 @@ target(0, 1) = -source(0, 1) * denominator; target(1, 0) = -source(1, 0) * denominator; target(1, 1) = source(0, 0) * denominator; + + return true; } else if (source.size1() == 3) { @@ -565,6 +568,8 @@ target(2, 0) = (d * h - e * g) * denominator; target(2, 1) = -(a * h - b * g) * denominator; target(2, 2) = (a * e - b * d) * denominator; + + return true; } else { @@ -581,12 +586,26 @@ if (boost::numeric::ublas::lu_factorize(a, permutation) != 0) { - LOG(ERROR) << "Cannot invert singular matrix"; - throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); + return false; + } + else + { + target = boost::numeric::ublas::identity_matrix(source.size1()); + lu_substitute(a, permutation, target); + return true; } + } + } + - target = boost::numeric::ublas::identity_matrix(source.size1()); - lu_substitute(a, permutation, target); + + void InvertMatrix(Matrix& target, + const Matrix& source) + { + if (!InvertMatrixUnsafe(target, source)) + { + LOG(ERROR) << "Cannot invert singular matrix"; + throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); } } diff -r 62964bfbbc00 -r 7105e51e4907 Framework/Toolbox/LinearAlgebra.h --- a/Framework/Toolbox/LinearAlgebra.h Fri Mar 02 11:21:12 2018 +0100 +++ b/Framework/Toolbox/LinearAlgebra.h Fri Mar 02 12:35:27 2018 +0100 @@ -222,6 +222,10 @@ void InvertMatrix(Matrix& target, const Matrix& source); + // This is the same as "InvertMatrix()", but without exception + bool InvertMatrixUnsafe(Matrix& target, + const Matrix& source); + void CreateSkewSymmetric(Matrix& s, const Vector& v);