|
|
| (2 intermediate revisions by one other user not shown) |
| Line 1: |
Line 1: |
| ITK optimizers are generic classes, which can be used independently of registration. This example demonstrates use of the itk::LevenbergMarquardtOptimizer class in optimizing a curve corrupted by Gaussian noise. | | {{warning|1=The media wiki content on this page is no longer maintained. The examples presented on the https://itk.org/Wiki/* pages likely require ITK version 4.13 or earlier releases. In many cases, the examples on this page no longer conform to the best practices for modern ITK versions. |
| | }} |
|
| |
|
| Contributed by: Davis Vigneault
| | [https://itk.org/ITKExamples[ITK Sphinx Examples]] |
| | |
| ==LMOptimization.cxx==
| |
| <source lang="cpp">
| |
| // Include the Levenberg-Marquardt optimizer and a custom cost function
| |
| #include "itkLevenbergMarquardtOptimizer.h"
| |
| #include "itkExampleCostFunction.h"
| |
| | |
| // Typedef the optimizer and cost function, for convenience | |
| typedef itk::LevenbergMarquardtOptimizer OptimizerType;
| |
| typedef itk::ExampleCostFunction CostType;
| |
| | |
| int main(int argc, char ** argv)
| |
| {
| |
| | |
| // Instantiate the cost function and optimizer
| |
| CostType::Pointer cost = CostType::New();
| |
| OptimizerType::Pointer optimizer = OptimizerType::New();
| |
| | |
| optimizer->SetNumberOfIterations( 100 );
| |
| optimizer->UseCostFunctionGradientOff();
| |
| optimizer->SetCostFunction( cost );
| |
| | |
| // This is the initial guess for the parameter values, which we set to one
| |
| CostType::ParametersType p(cost->GetNumberOfParameters());
| |
| p.Fill( 1 );
| |
| optimizer->SetInitialPosition(p);
| |
| optimizer->StartOptimization();
| |
| | |
| // Print out some information about the optimization
| |
| // The parameters come out to be near to, but not exactly [5.5, 0.5]
| |
| std::cout << "Position: " << optimizer->GetCurrentPosition() << std::endl;
| |
| std::cout << "Value: " << optimizer->GetValue() << std::endl;
| |
|
| |
| return EXIT_SUCCESS;
| |
| | |
| }
| |
| | |
| </source>
| |
| | |
| ==itkExampleCostFunction.h==
| |
| <source lang="cpp">
| |
| #ifndef itkExampleCostFunction_h
| |
| #define itkExampleCostFunction_h
| |
| | |
| #include <random>
| |
| #include <cmath>
| |
| #include <vector>
| |
| #include "itkMultipleValuedCostFunction.h"
| |
| | |
| namespace itk
| |
| {
| |
| class ExampleCostFunction :
| |
| public MultipleValuedCostFunction
| |
| {
| |
| public:
| |
| /** Standard class typedefs. */
| |
| typedef ExampleCostFunction Self;
| |
| typedef MultipleValuedCostFunction Superclass;
| |
| typedef SmartPointer<Self> Pointer;
| |
| typedef SmartPointer<const Self> ConstPointer;
| |
| | |
| /** Method for creation through the object factory. */
| |
| itkNewMacro(Self);
| |
| | |
| /** Run-time type information (and related methods). */
| |
| itkTypeMacro(ExampleCostFunction, MultipleValuedCostFunction);
| |
| | |
| // The equation we're fitting is y=C*e^(K*x)
| |
| // The free parameters which we're trying to fit are C and K
| |
| // Therefore, there are two parameters
| |
| unsigned int GetNumberOfParameters() const { return 2; }
| |
|
| |
| // We will take a curve with concrete values for C and K,
| |
| // which has been corrupted by Gaussian noise, and sample
| |
| // it at 100 points on the interval [0,1]. Each of these
| |
| // points will produce a residual with the expected value.
| |
| // Therefore, there are 100 values (aka residuals).
| |
| unsigned int GetNumberOfValues() const { return 100; }
| |
| | |
| // Calculate the residual array, given a set of parameters.
| |
| // We take parameters[0] to be C and parameters[1] to be K.
| |
| // Therefore, this is a matter of calculating the value of y
| |
| // at each of the sampled points, given the provided guesses
| |
| // for C and K, and returning the difference from the data.
| |
| MeasureType GetValue(const ParametersType ¶meters) const
| |
| {
| |
| MeasureType residual(this->GetNumberOfValues());
| |
| double predictedC = parameters[0];
| |
| double predictedK = parameters[1];
| |
| for (unsigned int i = 0; i < 100; ++i)
| |
| {
| |
| double position = double(i)/100;
| |
| double prediction = predictedC*exp(position*predictedK);
| |
| residual[i] = prediction - y[i];
| |
| }
| |
| return residual;
| |
| }
| |
| | |
| // The "derivative" is the Jacobian, which takes the derivative
| |
| // of each residual with respect to each parameter. Since this
| |
| // class does not provide a derivative, any optimizer using this
| |
| // cost function must be told explicitly not to ask for derivative,
| |
| // otherwise an exception will the thrown.
| |
| void GetDerivative(const ParametersType ¶meters, DerivativeType & derivative ) const
| |
| {
| |
| throw ExceptionObject(__FILE__,__LINE__,"No derivative available.");
| |
| }
| |
| | |
| protected:
| |
| ExampleCostFunction()
| |
| {
| |
| // Create some data
| |
| // Take the curve y = C*e^(K*x), add noise, and sample at 100 points on [0,1]
| |
| // C and K are our parameters
| |
| // In the actual data, these parameters are 5.5 and 0.5
| |
| typename std::default_random_engine generator;
| |
| typename std::normal_distribution<double> distribution(0.0,0.5);
| |
| double actualC = 5.5;
| |
| double actualK = 0.5;
| |
| for (unsigned int i = 0; i < 100; ++i)
| |
| {
| |
| double position = double(i)/100;
| |
| this->x.push_back(position);
| |
| this->y.push_back(actualC*exp(position*actualK)+distribution(generator));
| |
| }
| |
| };
| |
| ~ExampleCostFunction(){};
| |
| | |
| private:
| |
| ExampleCostFunction(const Self &); //purposely not implemented
| |
| void operator = (const Self &); //purposely not implemented
| |
| | |
| // The x and y positions of the data, created in the constructor
| |
| typename std::vector<double> x;
| |
| typename std::vector<double> y;
| |
| | |
| };
| |
| | |
| } // end namespace itk
| |
| | |
| #endif
| |
| </source>
| |
| {{ITKCMakeLists|{{SUBPAGENAME}}}}
| |
