Mercurial > hg > orthanc-stone
comparison Framework/Toolbox/LinearAlgebra.cpp @ 163:8c5b24892ed2 wasm
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 } |