@@ -25,23 +25,28 @@ LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
2525FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
2626DEALINGS IN THE SOFTWARE.
2727
28- Author : Moritz Dannhauer
29- Last modification : March 16 2014
30- ToDo: Padding is always enabled because of exit() in cleaver lib
28+ Author : Jaume Coll-Font
29+ Last modification : April 20 2017
3130*/
31+
3232#include < boost/bind.hpp>
3333#include < boost/lexical_cast.hpp>
3434
35+ // Tikhonov specific headers
3536#include < Core/Algorithms/Legacy/Inverse/TikhonovAlgoAbstractBase.h>
37+ #include < Core/Algorithms/Legacy/Inverse/TikhonovImpl.h>
38+ #include < Core/Algorithms/Legacy/Inverse/SolveInverseProblemWithStandardTikhonovImpl.h>
39+ #include < Core/Algorithms/Legacy/Inverse/SolveInverseProblemWithTikhonovSVD_impl.h>
3640
41+ // Datatypes
3742#include < Core/Datatypes/Matrix.h>
3843#include < Core/Datatypes/DenseMatrix.h>
3944#include < Core/Datatypes/DenseColumnMatrix.h>
4045#include < Core/Datatypes/SparseRowMatrix.h>
4146#include < Core/Datatypes/MatrixTypeConversions.h>
4247
48+ // SCIRun structural
4349#include < Core/Algorithms/Base/AlgorithmPreconditions.h>
44-
4550#include < Core/Logging/LoggerInterface.h>
4651#include < Core/Utils/Exception.h>
4752
@@ -61,7 +66,7 @@ const AlgorithmOutputName TikhonovAlgoAbstractBase::InverseSolution("InverseSolu
6166const AlgorithmOutputName TikhonovAlgoAbstractBase::RegularizationParameter (" RegularizationParameter" );
6267const AlgorithmOutputName TikhonovAlgoAbstractBase::RegInverse (" RegInverse" );
6368
64-
69+ const AlgorithmParameterName TikhonovAlgoAbstractBase::TikhonovImplementation ( " NoMethodSelected " );
6570const AlgorithmParameterName TikhonovAlgoAbstractBase::RegularizationMethod (" lambdaMethodComboBox" );
6671const AlgorithmParameterName TikhonovAlgoAbstractBase::regularizationChoice (" autoRadioButton" );
6772const AlgorithmParameterName TikhonovAlgoAbstractBase::LambdaFromDirectEntry (" lambdaDoubleSpinBox" );
@@ -109,7 +114,6 @@ TikhonovAlgoAbstractBase::TikhonovAlgoAbstractBase()
109114// //// CHECK IF INPUT MATRICES HAVE THE CORRECT SIZE
110115bool TikhonovAlgoAbstractBase::checkInputMatrixSizes ( const AlgorithmInput & input) const
111116{
112-
113117 // get inputs
114118 auto forwardMatrix_ = input.get <Matrix>(ForwardMatrix);
115119 auto measuredData_ = input.get <Matrix>(MeasuredPotentials);
@@ -200,16 +204,37 @@ AlgorithmOutput TikhonovAlgoAbstractBase::run(const AlgorithmInput & input) cons
200204
201205 // get Parameters
202206 auto RegularizationMethod_gotten = get (RegularizationMethod).toString ();
203-
204- // preAlocateInverseMatrices(forwardMatrix_,measuredData_,sourceWeighting_,sensorWeighting_);
205-
206- const int M = forwardMatrix_->nrows ();
207+ auto TikhonovImplementation_gotten = get (TikhonovImplementation).toString ();
207208
208209 // Alocate Variable
209210 DenseMatrix solution;
210211 double lambda_sq = 0 ;
211212 double lambda_ = 0 ;
212213
214+ // check input MATRICES
215+ checkInputMatrixSizes ( input );
216+
217+ // Determine specific Tikhonov Implementation
218+ TikhonovImpl *algoImpl;
219+ if ( TikhonovImplementation_gotten == " standardTikhonov" ){
220+ // get Parameters
221+ int regularizationChoice_ = get (regularizationChoice).toInt ();
222+ int regularizationSolutionSubcase_ = get (regularizationSolutionSubcase).toInt ();
223+ int regularizationResidualSubcase_ = get (regularizationResidualSubcase).toInt ();
224+
225+ algoImpl = new SolveInverseProblemWithStandardTikhonovImpl ( *castMatrix::toDense (forwardMatrix_), *castMatrix::toDense (measuredData_), *castMatrix::toDense (sourceWeighting_), *castMatrix::toDense (sensorWeighting_), regularizationChoice_, regularizationSolutionSubcase_, regularizationResidualSubcase_);
226+ }
227+ else if ( TikhonovImplementation_gotten == " TikhonovSVD" ){
228+ algoImpl = new SolveInverseProblemWithTikhonovSVD_impl ( *castMatrix::toDense (forwardMatrix_), *castMatrix::toDense (measuredData_), *castMatrix::toDense (sourceWeighting_), *castMatrix::toDense (sensorWeighting_));
229+ }
230+ else if ( TikhonovImplementation_gotten == " TikhonovTSVD" ){
231+ THROW_ALGORITHM_PROCESSING_ERROR (" Tikhonov TSVD not implemented yet" );
232+ }
233+ else {
234+ THROW_ALGORITHM_PROCESSING_ERROR (" Not a valid Tikhonov Implementation selection" );
235+ }
236+
237+ // preAlocateInverseMatrices(forwardMatrix_,measuredData_,sourceWeighting_,sensorWeighting_);
213238
214239 // Get Regularization parameter(s) : Lambda
215240 if ((RegularizationMethod_gotten == " single" ) || (RegularizationMethod_gotten == " slider" ))
@@ -227,7 +252,7 @@ AlgorithmOutput TikhonovAlgoAbstractBase::run(const AlgorithmInput & input) cons
227252 }
228253 else if (RegularizationMethod_gotten == " lcurve" )
229254 {
230- lambda_ = computeLcurve ( input );
255+ lambda_ = computeLcurve ( algoImpl, input );
231256 }
232257 else
233258 {
@@ -239,7 +264,7 @@ AlgorithmOutput TikhonovAlgoAbstractBase::run(const AlgorithmInput & input) cons
239264
240265
241266 // compute inverse solution
242- solution = computeInverseSolution (lambda_sq, true );
267+ solution = algoImpl-> computeInverseSolution (lambda_sq, true );
243268
244269 //
245270 // // set final result
@@ -266,7 +291,7 @@ AlgorithmOutput TikhonovAlgoAbstractBase::run(const AlgorithmInput & input) cons
266291
267292// /////////////////////////
268293// ///// compute L-curve
269- double TikhonovAlgoAbstractBase::computeLcurve (const AlgorithmInput & input ) const
294+ double TikhonovAlgoAbstractBase::computeLcurve ( const SCIRun::Core::Algorithms::Inverse::TikhonovImpl * algoImpl, const AlgorithmInput & input ) const
270295{
271296
272297 // get inputs
@@ -309,8 +334,8 @@ double TikhonovAlgoAbstractBase::computeLcurve(const AlgorithmInput & input ) co
309334 // set current lambda
310335 lambda_sq = lambdaArray[j] * lambdaArray[j];
311336
312- // COMPUTE INVERSE SOLUTION // Todo: @JCOLLFONT function needs to be defined
313- solution = computeInverseSolution ( lambda_sq, false );
337+ // COMPUTE INVERSE SOLUTION
338+ solution = algoImpl-> computeInverseSolution ( lambda_sq, false );
314339
315340
316341 // if using source regularization matrix, apply it to compute Rx (for the eta computations)
0 commit comments