Skip to content

Commit 2507820

Browse files
jcollfontjcollfont
committed
created the basic modules for the Tikhonov TSVD module. It does not work properly. The assumption that the possible lambdas are integers
1 parent 97e7aad commit 2507820

12 files changed

+910
-2
lines changed

‎src/Core/Algorithms/Legacy/Inverse/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,13 +30,15 @@ SET(Algorithms_Legacy_Inverse_SRCS
3030
TikhonovAlgoAbstractBase.cc
3131
SolveInverseProblemWithStandardTikhonovImpl.cc
3232
SolveInverseProblemWithTikhonovSVD_impl.cc
33+
SolveInverseProblemWithTikhonovTSVD_impl.cc
3334
)
3435

3536
SET(Algorithms_Legacy_Inverse_HEADERS
3637
TikhonovAlgoAbstractBase.h
3738
TikhonovImpl.h
3839
SolveInverseProblemWithStandardTikhonovImpl.h
3940
SolveInverseProblemWithTikhonovSVD_impl.h
41+
SolveInverseProblemWithTikhonovTSVD_impl.h
4042
share.h
4143
)
4244

‎src/Core/Algorithms/Legacy/Inverse/SolveInverseProblemWithTikhonovTSVD_impl.cc‎

Lines changed: 143 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,143 @@
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+
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_)
93+
{
94+
95+
// Compute the SVD of the forward matrix
96+
Eigen::JacobiSVD<SCIRun::Core::Datatypes::DenseMatrix::EigenBase> SVDdecomposition( forwardMatrix_, Eigen::ComputeFullU | Eigen::ComputeFullV);
97+
98+
// alocate the left and right singular vectors and the singular values
99+
svd_MatrixU = SVDdecomposition.matrixU();
100+
svd_MatrixV = SVDdecomposition.matrixV();
101+
svd_SingularValues = SVDdecomposition.singularValues();
102+
103+
// determine rank
104+
rank = SVDdecomposition.nonzeroSingularValues();
105+
106+
// Compute the projection of data y on the left singular vectors
107+
Uy = svd_MatrixU.transpose() * (measuredData_);
108+
}
109+
110+
//////////////////////////////////////////////////////////////////////
111+
// THIS FUNCTION returns regularized solution by tikhonov method
112+
//////////////////////////////////////////////////////////////////////
113+
SCIRun::Core::Datatypes::DenseMatrix SolveInverseProblemWithTikhonovTSVD_impl::computeInverseSolution( double truncationPoint, bool inverseCalculation ) const
114+
{
115+
116+
// prealocate matrices
117+
const int N = svd_MatrixV.cols();
118+
const int M = svd_MatrixU.rows();
119+
const int numTimeSamples = Uy.ncols();
120+
DenseMatrix solution(DenseMatrix::Zero(N,numTimeSamples));
121+
DenseMatrix tempInverse(DenseMatrix::Zero(N,M));
122+
123+
// Compute inverse SolveInverseProblemWithTikhonovTSVD
124+
for (int rr=0; rr < truncationPoint ; rr++)
125+
{
126+
// evaluate filter factor
127+
double singVal = svd_SingularValues[rr];
128+
double filterFactor_i = 1 / ( singVal ) * Uy(rr);
129+
130+
// u[date solution
131+
solution += filterFactor_i * svd_MatrixV.col(rr);
132+
133+
// update inverse operator
134+
if (inverseCalculation)
135+
tempInverse += filterFactor_i * ( svd_MatrixV.col(rr) * svd_MatrixU.col(rr).transpose() );
136+
}
137+
138+
// output solutions
139+
// if (inverseCalculation)
140+
// inverseMatrix_.reset( new SCIRun::Core::Datatypes::DenseMatrix(tempInverse) );
141+
142+
return solution;
143+
}

‎src/Core/Algorithms/Legacy/Inverse/SolveInverseProblemWithTikhonovTSVD_impl.h‎

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
/*
2+
For more information, please see: http://software.sci.utah.edu
3+
4+
The MIT License
5+
6+
Copyright (c) 2015 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 : SolveInverseProblemWithTikhonov.cc
30+
// Author : Moritz Dannhauer and Ayla Khan
31+
// Date : 15 Aug. 2012
32+
33+
#ifndef BioPSE_SolveInverseProblemWithTikhonovTSVDimpl_H__
34+
#define BioPSE_SolveInverseProblemWithTikhonovTSVDimpl_H__
35+
36+
#include <vector>
37+
38+
#include <boost/utility.hpp>
39+
#include <boost/function.hpp>
40+
41+
#include <Core/Datatypes/MatrixFwd.h>
42+
#include <Core/Datatypes/DenseMatrix.h>
43+
#include <Core/Datatypes/DenseColumnMatrix.h>
44+
#include <Core/Logging/LoggerFwd.h>
45+
46+
#include <Core/Algorithms/Legacy/Inverse/TikhonovImpl.h>
47+
48+
#include <Core/Algorithms/Legacy/Inverse/share.h>
49+
50+
51+
namespace SCIRun {
52+
namespace Core {
53+
namespace Algorithms {
54+
namespace Inverse {
55+
56+
class SCISHARE SolveInverseProblemWithTikhonovTSVD_impl : public TikhonovImpl
57+
{
58+
59+
60+
public:
61+
SolveInverseProblemWithTikhonovTSVD_impl(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_)
62+
{
63+
preAlocateInverseMatrices( forwardMatrix_, measuredData_ , sourceWeighting_, sensorWeighting_, matrixU_, singularValues_, matrixV_ );
64+
};
65+
66+
SolveInverseProblemWithTikhonovTSVD_impl(const SCIRun::Core::Datatypes::DenseMatrix& forwardMatrix_, const SCIRun::Core::Datatypes::DenseMatrix& measuredData_ , const SCIRun::Core::Datatypes::DenseMatrix& sourceWeighting_, const SCIRun::Core::Datatypes::DenseMatrix& sensorWeighting_)
67+
{
68+
preAlocateInverseMatrices( forwardMatrix_, measuredData_ , sourceWeighting_, sensorWeighting_);
69+
};
70+
71+
private:
72+
73+
// Data Members
74+
int rank;
75+
SCIRun::Core::Datatypes::DenseMatrix svd_MatrixU;
76+
SCIRun::Core::Datatypes::DenseColumnMatrix svd_SingularValues;
77+
SCIRun::Core::Datatypes::DenseMatrix svd_MatrixV;
78+
79+
SCIRun::Core::Datatypes::DenseMatrix Uy;
80+
81+
// Methods
82+
void 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_);
83+
void 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_);
84+
85+
virtual SCIRun::Core::Datatypes::DenseMatrix computeInverseSolution( double truncationPoint, bool inverseCalculation) const;
86+
// bool checkInputMatrixSizes(); // DEFINED IN PARENT, MIGHT WANT TO OVERRIDE SOME OTHER TIME
87+
88+
89+
};
90+
}}}}
91+
92+
#endif

‎src/Core/Algorithms/Legacy/Inverse/TikhonovAlgoAbstractBase.cc‎

Lines changed: 18 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@ Last modification : April 20 2017
3737
#include <Core/Algorithms/Legacy/Inverse/TikhonovImpl.h>
3838
#include <Core/Algorithms/Legacy/Inverse/SolveInverseProblemWithStandardTikhonovImpl.h>
3939
#include <Core/Algorithms/Legacy/Inverse/SolveInverseProblemWithTikhonovSVD_impl.h>
40+
#include <Core/Algorithms/Legacy/Inverse/SolveInverseProblemWithTikhonovTSVD_impl.h>
4041

4142
// Datatypes
4243
#include <Core/Datatypes/Matrix.h>
@@ -236,15 +237,30 @@ AlgorithmOutput TikhonovAlgoAbstractBase::run(const AlgorithmInput & input) cons
236237
auto singularValues_ = input.get<Matrix>(TikhonovAlgoAbstractBase::singularValues);
237238
auto matrixV_ = input.get<Matrix>(TikhonovAlgoAbstractBase::matrixV);
238239

239-
// If there is a missing matrix from the precomputed SVD input
240+
// If there is a missing matrix from the precomputed SVD input
240241
if ( (matrixU_ == NULL) || (singularValues_ == NULL) || ( matrixV_ == NULL) )
241242
algoImpl = new SolveInverseProblemWithTikhonovSVD_impl( *castMatrix::toDense(forwardMatrix_), *castMatrix::toDense(measuredData_), *castMatrix::toDense(sourceWeighting_), *castMatrix::toDense(sensorWeighting_) );
242243
else
243244
algoImpl = new SolveInverseProblemWithTikhonovSVD_impl( *castMatrix::toDense(forwardMatrix_), *castMatrix::toDense(measuredData_), *castMatrix::toDense(sourceWeighting_), *castMatrix::toDense(sensorWeighting_), *castMatrix::toDense(matrixU_), *castMatrix::toDense(singularValues_), *castMatrix::toDense(matrixV_) );
244245

245246
}
246247
else if ( TikhonovImplementation_gotten == std::string("TikhonovTSVD") ){
247-
THROW_ALGORITHM_PROCESSING_ERROR("Tikhonov TSVD not implemented yet");
248+
249+
// get Parameters
250+
int regularizationChoice_ = get(Parameters::regularizationChoice).toInt();
251+
int regularizationSolutionSubcase_ = get(Parameters::regularizationSolutionSubcase).toInt();
252+
int regularizationResidualSubcase_ = get(Parameters::regularizationResidualSubcase).toInt();
253+
254+
// get TikhonovSVD special inputs
255+
auto matrixU_ = input.get<Matrix>(TikhonovAlgoAbstractBase::matrixU);
256+
auto singularValues_ = input.get<Matrix>(TikhonovAlgoAbstractBase::singularValues);
257+
auto matrixV_ = input.get<Matrix>(TikhonovAlgoAbstractBase::matrixV);
258+
259+
// If there is a missing matrix from the precomputed SVD input
260+
if ( (matrixU_ == NULL) || (singularValues_ == NULL) || ( matrixV_ == NULL) )
261+
algoImpl = new SolveInverseProblemWithTikhonovTSVD_impl( *castMatrix::toDense(forwardMatrix_), *castMatrix::toDense(measuredData_), *castMatrix::toDense(sourceWeighting_), *castMatrix::toDense(sensorWeighting_) );
262+
else
263+
algoImpl = new SolveInverseProblemWithTikhonovTSVD_impl( *castMatrix::toDense(forwardMatrix_), *castMatrix::toDense(measuredData_), *castMatrix::toDense(sourceWeighting_), *castMatrix::toDense(sensorWeighting_), *castMatrix::toDense(matrixU_), *castMatrix::toDense(singularValues_), *castMatrix::toDense(matrixV_) );
248264
}
249265
else{
250266
THROW_ALGORITHM_PROCESSING_ERROR("Not a valid Tikhonov Implementation selection");

‎src/Interface/Modules/Inverse/CMakeLists.txt‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,17 +30,20 @@
3030
SET(Interface_Modules_Inverse_FORMS
3131
SolveInverseProblemWithTikhonovSVDDialog.ui
3232
SolveInverseProblemWithTikhonov.ui
33+
SolveInverseProblemWithTikhonovTSVDDialog.ui
3334
)
3435

3536
SET(Interface_Modules_Inverse_HEADERS
3637
SolveInverseProblemWithTikhonovSVDDialog.h
3738
SolveInverseProblemWithTikhonovDialog.h
39+
SolveInverseProblemWithTikhonovTSVDDialog.h
3840
share.h
3941
)
4042

4143
SET(Interface_Modules_Inverse_SOURCES
4244
SolveInverseProblemWithTikhonovSVDDialog.cc
4345
SolveInverseProblemWithTikhonovDialog.cc
46+
SolveInverseProblemWithTikhonovTSVDDialog.cc
4447
)
4548

4649
QT4_WRAP_UI(Interface_Modules_Inverse_FORMS_HEADERS ${Interface_Modules_Inverse_FORMS})

0 commit comments

Comments
 (0)