Mercurial > hg > orthanc-stone

comparison Framework/Toolbox/LinearAlgebra.cpp @ 163:8c5b24892ed2 wasm

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
LinearAlgebra::InvertMatrix
author Sebastien Jodogne <s.jodogne@gmail.com>
date Wed, 14 Feb 2018 14:26:26 +0100
parents 197a5ddaf68c
children 432b1f812d14
comparison
equal deleted inserted replaced
162:77715c340767 163:8c5b24892ed2
447 !IsCloseToZero(r(2, 1))) 447 !IsCloseToZero(r(2, 1)))
448 { 448 {
449 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); 449 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
450 } 450 }
451 } 451 }
452
453
454 void InvertMatrix(Matrix& target,
455 const Matrix& source)
456 {
457 if (source.size1() != source.size2())
458 {
459 LOG(ERROR) << "Inverse only exists for square matrices";
460 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
461 }
462
463 double determinant = ComputeDeterminant(source);
464
465 if (IsCloseToZero(determinant))
466 {
467 LOG(ERROR) << "Cannot invert singular matrix";
468 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
469 }
470
471 double denominator = 1.0 / determinant;
472
473 target.resize(source.size1(), source.size2());
474
475 if (source.size1() == 1)
476 {
477 target(0, 0) = denominator;
478 }
479 else if (source.size1() == 2)
480 {
481 // https://en.wikipedia.org/wiki/Invertible_matrix#Inversion_of_2_%C3%97_2_matrices
482 target(0, 0) = source(1, 1) * denominator;
483 target(0, 1) = -source(0, 1) * denominator;
484 target(1, 0) = -source(1, 0) * denominator;
485 target(1, 1) = source(0, 0) * denominator;
486 }
487 else if (source.size1() == 3)
488 {
489 // https://en.wikipedia.org/wiki/Invertible_matrix#Inversion_of_3_%C3%97_3_matrices
490 const double a = source(0, 0);
491 const double b = source(0, 1);
492 const double c = source(0, 2);
493 const double d = source(1, 0);
494 const double e = source(1, 1);
495 const double f = source(1, 2);
496 const double g = source(2, 0);
497 const double h = source(2, 1);
498 const double i = source(2, 2);
499
500 target(0, 0) = (e * i - f * h) * denominator;
501 target(0, 1) = -(b * i - c * h) * denominator;
502 target(0, 2) = (b * f - c * e) * denominator;
503 target(1, 0) = -(d * i - f * g) * denominator;
504 target(1, 1) = (a * i - c * g) * denominator;
505 target(1, 2) = -(a * f - c * d) * denominator;
506 target(2, 0) = (d * h - e * g) * denominator;
507 target(2, 1) = -(a * h - b * g) * denominator;
508 target(2, 2) = (a * e - b * d) * denominator;
509 }
510 else
511 {
512 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
513 }
514 }
452 } 515 }
453 } 516 }