added functionality to allow each inverse method to generate its special sequence of lambdas. If not implemented, the default option is logarithmic increase

Author: jcollfont

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

@@ -144,3 +144,19 @@ SCIRun::Core::Datatypes::DenseMatrix SolveInverseProblemWithTikhonovTSVD_impl::c
 
         return solution;
 }
+
+//////////////////////////////////////////////////////////////////////
+// THIS FUNCTION returns a string of lambdas from which the L-curve is computed
+//////////////////////////////////////////////////////////////////////
+std::vector<double> SolveInverseProblemWithTikhonovTSVD_impl::computeLambdaArray( double lambdaMin, double lambdaMax, int nLambda ) const
+{
+	std::vector<double> lambdaArray(nLambda,0.0);
+	const double lam_step = 1;
+
+	lambdaArray[0] = lambdaMin;
+	for (int j = 1; j < nLambda; j++)
+	{
+		lambdaArray[j] = lambdaArray[j-1]  + lam_step;
+	}
+	return lambdaArray;
+}

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

@@ -83,6 +83,7 @@ namespace Inverse {
 				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_);
 
 		        virtual SCIRun::Core::Datatypes::DenseMatrix computeInverseSolution( double truncationPoint, bool inverseCalculation) const;
+				std::vector<double> computeLambdaArray( double lambdaMin, double lambdaMax, int nLambda ) const;
 		        //      bool checkInputMatrixSizes(); // DEFINED IN PARENT, MIGHT WANT TO OVERRIDE SOME OTHER TIME
 
 

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

@@ -336,15 +336,17 @@ double TikhonovAlgoAbstractBase::computeLcurve( const SCIRun::Core::Algorithms::
     const int nLambda = get(Parameters::LambdaNum).toInt();
 	const double lambdaMin_ = get(Parameters::LambdaMin).toDouble();
 	const double lambdaMax_ = get(Parameters::LambdaMax).toDouble();
-    const double lam_step = (log10(lambdaMax_) - log10(lambdaMin_))  / (nLambda-1);
-	std::cout << "Lambda power step: " << lam_step << ". Number: "<< nLambda <<". Lambda min: " << lambdaMin_ << ". Lambda max: "<< lambdaMax_<< ". Ratio: "<<  lambdaMax_ / lambdaMin_ << std::endl;
-    double lambda = 0;
+	double lambda = 0;
 
-    // prealocate vector of lambdas and eta and rho
-    std::vector<double> lambdaArray(nLambda, 0.0);
+	// prealocate vector of lambdas and eta and rho
     std::vector<double> rho(nLambda, 0.0);
     std::vector<double> eta(nLambda, 0.0);
 
+	// Compute the lambda array
+	std::vector<double> lambdaArray(nLambda, 0.0);
+    lambdaArray = algoImpl->computeLambdaArray( lambdaMin_, lambdaMax_, nLambda );
+
+
     DenseMatrix CAx, Rx;
     DenseMatrix solution;
 
@@ -353,10 +355,6 @@ double TikhonovAlgoAbstractBase::computeLcurve( const SCIRun::Core::Algorithms::
     // for all lambdas
     for (int j = 0; j < nLambda; j++)
     {
-        if (j)
-        {
-            lambdaArray[j] = lambdaArray[j-1] * pow(10.0,lam_step);
-        }
 
         // COMPUTE INVERSE SOLUTION
         solution = algoImpl->computeInverseSolution( lambdaArray[j], false);
@@ -386,7 +384,7 @@ double TikhonovAlgoAbstractBase::computeLcurve( const SCIRun::Core::Algorithms::
             CAx = residualSolution;
 
 
-        // compute rho and eta
+        // compute rho and eta. Using Frobenious norm when using matrices
 		rho[j] = CAx.norm();
 		eta[j] = Rx.norm();
 

src/Core/Algorithms/Legacy/Inverse/TikhonovImpl.h

@@ -48,6 +48,21 @@ namespace Inverse {
 
 		virtual SCIRun::Core::Datatypes::DenseMatrix computeInverseSolution( double lambda_sq, bool inverseCalculation) const = 0;
 
+		// default lambda step. Can ve overriden if necessary (see TSVD as reference)
+		virtual std::vector<double> computeLambdaArray( double lambdaMin, double lambdaMax, int nLambda ) const {
+
+			std::vector<double> lambdaArray(nLambda,0.0);
+			double lam_step = (log10(lambdaMax) - log10(lambdaMin)) / (nLambda-1);
+
+			lambdaArray[0] = lambdaMin;
+			for (int j = 1; j < nLambda; j++)
+		    {
+		        lambdaArray[j] = lambdaArray[j-1] * pow(10.0,lam_step);
+			}
+
+			return lambdaArray;
+		};
+
 	};
 
 	}}}}

src/Modules/Legacy/Inverse/SolveInverseProblemWithTikhonovTSVD.cc

@@ -122,7 +122,7 @@ void SolveInverseProblemWithTikhonovTSVD::execute()
 		// set parameters
 		auto state = get_state();
 
-		// state->setValue( Parameters::TikhonovImplementation, std::string("TikhonovTSVD") );
+		state->setValue( Parameters::TikhonovImplementation, std::string("TikhonovTSVD") );
 		setAlgoStringFromState(Parameters::TikhonovImplementation);
 		setAlgoOptionFromState(Parameters::RegularizationMethod);
 		setAlgoDoubleFromState(Parameters::LambdaFromDirectEntry);