Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkBiExpFitFunctor.cpp
Go to the documentation of this file.
1 /*===================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
7 All rights reserved.
8 
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
11 A PARTICULAR PURPOSE.
12 
13 See LICENSE.txt or http://www.mitk.org for details.
14 
15 ===================================================================*/
16 
17 #include "itkBiExpFitFunctor.h"
18 #include <math.h>
19 #include <iostream>
20 #include <iomanip>
21 
22 void itk::BiExpFitFunctor::operator()(vnl_matrix<double> & newSignal,const vnl_matrix<double> & SignalMatrix, const double & S0)
23 {
24 
25  vnl_vector<double> initalGuess(3);
26  // initialize Least Squres Function
27  // SignalMatrix.cols() defines the number of shells points
28  lestSquaresFunction model(SignalMatrix.cols());
29  model.set_bvalues(m_BValueList);// set BValue Vector e.g.: [1000, 2000, 3000] <- shell b Values
30 
31  // initialize Levenberg Marquardt
32  vnl_levenberg_marquardt minimizer(model);
33  minimizer.set_max_function_evals(1000); // Iterations
34  minimizer.set_f_tolerance(1e-10); // Function tolerance
35 
36  // for each Direction calculate LSF Coeffs ADC & AKC
37  for(unsigned int i = 0 ; i < SignalMatrix.rows(); i++)
38  {
39  model.set_measurements(SignalMatrix.get_row(i));
40  model.set_reference_measurement(S0);
41 
42  initalGuess.put(0, 0.f); // ADC_slow
43  initalGuess.put(1, 0.009f); // ADC_fast
44  initalGuess.put(2, 0.7f); // lambda
45 
46  // start Levenberg-Marquardt
47  minimizer.minimize_without_gradient(initalGuess);
48 
49  const double & ADC_slow = initalGuess.get(0);
50  const double & ADC_fast = initalGuess.get(1);
51  const double & lambda = initalGuess(2);
52 
53  newSignal.put(i, 0, S0 * (lambda * std::exp(-m_TargetBvalue * ADC_slow) + (1-lambda)* std::exp(-m_TargetBvalue * ADC_fast)));
54  newSignal.put(i, 1, minimizer.get_end_error()); // RMS Error
55 
56  //OUTPUT FOR EVALUATION
57  /*std::cout << std::scientific << std::setprecision(5)
58  << ADC_slow << "," // lambda
59  << ADC_fast << "," // alpha
60  << lambda << "," // lambda
61  << S0 << "," // S0 value
62  << minimizer.get_end_error() << ","; // End error
63  for(unsigned int j = 0; j < SignalMatrix.get_row(i).size(); j++ ){
64  std::cout << std::scientific << std::setprecision(5) << SignalMatrix.get_row(i)[j]; // S_n Values corresponding to shell 1 to shell n
65  if(j != SignalMatrix.get_row(i).size()-1) std::cout << ",";
66  }
67  std::cout << std::endl;*/
68  }
69 
70 }
The lestSquaresFunction struct for Non-Linear-Least-Squres fit of Biexponential model.
vnl_vector< double > m_BValueList
void set_bvalues(const vnl_vector< double > &x)
void operator()(vnl_matrix< double > &newSignal, const vnl_matrix< double > &SignalMatrix, const double &S0) override
operator ()