|
| 1 | +/* |
| 2 | + For more information, please see: http://software.sci.utah.edu |
| 3 | +
|
| 4 | + The MIT License |
| 5 | +
|
| 6 | + Copyright (c) 2009 Scientific Computing and Imaging Institute, |
| 7 | + University of Utah. |
| 8 | +
|
| 9 | +
|
| 10 | + Permission is hereby granted, free of charge, to any person obtaining a |
| 11 | + copy of this software and associated documentation files (the "Software"), |
| 12 | + to deal in the Software without restriction, including without limitation |
| 13 | + the rights to use, copy, modify, merge, publish, distribute, sublicense, |
| 14 | + and/or sell copies of the Software, and to permit persons to whom the |
| 15 | + Software is furnished to do so, subject to the following conditions: |
| 16 | +
|
| 17 | + The above copyright notice and this permission notice shall be included |
| 18 | + in all copies or substantial portions of the Software. |
| 19 | +
|
| 20 | + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS |
| 21 | + OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |
| 22 | + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL |
| 23 | + THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER |
| 24 | + LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING |
| 25 | + FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER |
| 26 | + DEALINGS IN THE SOFTWARE. |
| 27 | + */ |
| 28 | + |
| 29 | +// File : SolveInverseProblemWithTikhonovTSVD.cc |
| 30 | +// Author : Yesim Serinagaoglu & Alireza Ghodrati |
| 31 | +// Date : 07 Aug. 2001 |
| 32 | +// Last update: Dec 2011 |
| 33 | + |
| 34 | + |
| 35 | +// SCIRUN lybraries |
| 36 | +#include <boost/bind.hpp> |
| 37 | +#include <boost/lexical_cast.hpp> |
| 38 | +#include <Core/Datatypes/Matrix.h> |
| 39 | +#include <Core/Datatypes/DenseMatrix.h> |
| 40 | +#include <Core/Datatypes/DenseColumnMatrix.h> |
| 41 | +#include <Core/Datatypes/SparseRowMatrix.h> |
| 42 | +#include <Core/Datatypes/MatrixTypeConversions.h> |
| 43 | +#include <Core/Algorithms/Base/AlgorithmPreconditions.h> |
| 44 | +#include <Core/Logging/LoggerInterface.h> |
| 45 | +#include <Core/Utils/Exception.h> |
| 46 | + |
| 47 | +// Tikhonov inverse libraries |
| 48 | +#include <Core/Algorithms/Legacy/Inverse/TikhonovAlgoAbstractBase.h> |
| 49 | +#include <Core/Algorithms/Legacy/Inverse/SolveInverseProblemWithTikhonovTSVD_impl.h> |
| 50 | + |
| 51 | +// EIGEN LIBRARY |
| 52 | +#include <Eigen/Eigen> |
| 53 | +#include <Eigen/SVD> |
| 54 | + |
| 55 | + |
| 56 | +using namespace SCIRun; |
| 57 | +using namespace SCIRun::Core::Datatypes; |
| 58 | +// using namespace SCIRun::Modules::Inverse; |
| 59 | +// using namespace SCIRun::Dataflow::Networks; |
| 60 | +using namespace SCIRun::Core::Logging; |
| 61 | +using namespace SCIRun::Core::Algorithms; |
| 62 | +using namespace SCIRun::Core::Algorithms::Inverse; |
| 63 | + |
| 64 | + |
| 65 | + |
| 66 | +/////////////////////////////////////////////////////////////////// |
| 67 | +/////// prealocate Matrices for inverse compuation |
| 68 | +/// This function precalcualtes the SVD of the forward matrix and prepares singular vectors and values for posterior computations |
| 69 | +/////////////////////////////////////////////////////////////////// |
| 70 | +void SolveInverseProblemWithTikhonovTSVD_impl::preAlocateInverseMatrices(const SCIRun::Core::Datatypes::DenseMatrix& forwardMatrix_, const SCIRun::Core::Datatypes::DenseMatrix& measuredData_ , const SCIRun::Core::Datatypes::DenseMatrix& sourceWeighting_, const SCIRun::Core::Datatypes::DenseMatrix& sensorWeighting_, const SCIRun::Core::Datatypes::DenseMatrix& matrixU_, const SCIRun::Core::Datatypes::DenseMatrix& singularValues_, const SCIRun::Core::Datatypes::DenseMatrix& matrixV_) |
| 71 | +{ |
| 72 | + |
| 73 | + // alocate U and V matrices |
| 74 | + svd_MatrixU = matrixU_; |
| 75 | + svd_MatrixV = matrixV_; |
| 76 | + |
| 77 | + // alocate singular values |
| 78 | + if (singularValues_.ncols() == 1 ){ |
| 79 | + svd_SingularValues = singularValues_; |
| 80 | + } |
| 81 | + else{ |
| 82 | + svd_SingularValues = singularValues_.diagonal(); |
| 83 | + } |
| 84 | + |
| 85 | + // Compute the projection of data y on the left singular vectors |
| 86 | + Uy = svd_MatrixU.transpose() * (measuredData_); |
| 87 | + |
| 88 | + // determine rank |
| 89 | + rank = svd_SingularValues.nrows(); |
| 90 | + |
| 91 | +} |
| 92 | + |
| 93 | +void SolveInverseProblemWithTikhonovTSVD_impl::preAlocateInverseMatrices(const SCIRun::Core::Datatypes::DenseMatrix& forwardMatrix_, const SCIRun::Core::Datatypes::DenseMatrix& measuredData_ , const SCIRun::Core::Datatypes::DenseMatrix& sourceWeighting_, const SCIRun::Core::Datatypes::DenseMatrix& sensorWeighting_) |
| 94 | +{ |
| 95 | + |
| 96 | + // Compute the SVD of the forward matrix |
| 97 | + Eigen::JacobiSVD<SCIRun::Core::Datatypes::DenseMatrix::EigenBase> SVDdecomposition( forwardMatrix_, Eigen::ComputeFullU | Eigen::ComputeFullV); |
| 98 | + |
| 99 | + // alocate the left and right singular vectors and the singular values |
| 100 | + svd_MatrixU = SVDdecomposition.matrixU(); |
| 101 | + svd_MatrixV = SVDdecomposition.matrixV(); |
| 102 | + svd_SingularValues = SVDdecomposition.singularValues(); |
| 103 | + |
| 104 | + // determine rank |
| 105 | + rank = SVDdecomposition.nonzeroSingularValues(); |
| 106 | + |
| 107 | + // Compute the projection of data y on the left singular vectors |
| 108 | + Uy = svd_MatrixU.transpose() * (measuredData_); |
| 109 | +} |
| 110 | + |
| 111 | +////////////////////////////////////////////////////////////////////// |
| 112 | +// THIS FUNCTION returns regularized solution by tikhonov method |
| 113 | +////////////////////////////////////////////////////////////////////// |
| 114 | +SCIRun::Core::Datatypes::DenseMatrix SolveInverseProblemWithTikhonovTSVD_impl::computeInverseSolution( double lambda, bool inverseCalculation ) const |
| 115 | +{ |
| 116 | + |
| 117 | + // prealocate matrices |
| 118 | + const int N = svd_MatrixV.cols(); |
| 119 | + const int M = svd_MatrixU.rows(); |
| 120 | + const int numTimeSamples = Uy.ncols(); |
| 121 | + DenseMatrix solution(DenseMatrix::Zero(N,numTimeSamples)); |
| 122 | + DenseMatrix tempInverse(DenseMatrix::Zero(N,M)); |
| 123 | + |
| 124 | + const int truncationPoint = Min( int(lambda), rank, int(9999999999999) ); |
| 125 | + |
| 126 | + // Compute inverse SolveInverseProblemWithTikhonovTSVD |
| 127 | + for (int rr=0; rr < truncationPoint ; rr++) |
| 128 | + { |
| 129 | + // evaluate filter factor |
| 130 | + double singVal = svd_SingularValues[rr]; |
| 131 | + double filterFactor_i = 1 / ( singVal ) * Uy(rr); |
| 132 | + |
| 133 | + // u[date solution |
| 134 | + solution += filterFactor_i * svd_MatrixV.col(rr); |
| 135 | + |
| 136 | + // update inverse operator |
| 137 | + if (inverseCalculation) |
| 138 | + tempInverse += filterFactor_i * ( svd_MatrixV.col(rr) * svd_MatrixU.col(rr).transpose() ); |
| 139 | + } |
| 140 | + |
| 141 | + // output solutions |
| 142 | + // if (inverseCalculation) |
| 143 | + // inverseMatrix_.reset( new SCIRun::Core::Datatypes::DenseMatrix(tempInverse) ); |
| 144 | + |
| 145 | + return solution; |
| 146 | +} |
0 commit comments