ITK/Examples/Math/LMOptimization

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.

Contributed by: Davis Vigneault

LMOptimization.cxx

<source lang="cpp"> // Include the Levenberg-Marquardt optimizer and a custom cost function

1. include "itkLevenbergMarquardtOptimizer.h"
2. include "itkExampleCostFunction.h"

// Typedef the optimizer and cost function, for convenience typedef itk::LevenbergMarquardtOptimizer OptimizerType; typedef itk::ExampleCostFunction CostType;

int main(int, char *[]) {

 // 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">

1. ifndef itkExampleCostFunction_h
2. define itkExampleCostFunction_h

1. include "itkMultipleValuedCostFunction.h"
2. include "itkMersenneTwisterRandomVariateGenerator.h"
3. include <cmath>
4. include <vector>

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 &parameters) 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 &parameters, 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
   typedef itk::Statistics::MersenneTwisterRandomVariateGenerator
     GeneratorType;
   GeneratorType::Pointer randomEngine = GeneratorType::New();
   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) +
                       randomEngine->GetUniformVariate(0.0, 0.5));
     }
   };
 ~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
 std::vector<double> x;
 std::vector<double> y;

};

} // end namespace itk

1. endif

</source>

CMakeLists.txt

<syntaxhighlight lang="cmake"> cmake_minimum_required(VERSION 3.9.5)

project(LMOptimization)

find_package(ITK REQUIRED) include(${ITK_USE_FILE}) if (ITKVtkGlue_LOADED)

 find_package(VTK REQUIRED)
 include(${VTK_USE_FILE})

endif()

add_executable(LMOptimization MACOSX_BUNDLE LMOptimization.cxx)

if( "${ITK_VERSION_MAJOR}" LESS 4 )

 target_link_libraries(LMOptimization ITKReview ${ITK_LIBRARIES})

else( "${ITK_VERSION_MAJOR}" LESS 4 )

 target_link_libraries(LMOptimization ${ITK_LIBRARIES})

endif( "${ITK_VERSION_MAJOR}" LESS 4 )

</syntaxhighlight>

Download and Build LMOptimization

Click here to download LMOptimization and its CMakeLists.txt file. Once the tarball LMOptimization.tar has been downloaded and extracted,

cd LMOptimization/build

• If ITK is installed:

cmake ..

• If ITK is not installed but compiled on your system, you will need to specify the path to your ITK build:

cmake -DITK_DIR:PATH=/home/me/itk_build ..

Build the project:

make

and run it:

./LMOptimization

WINDOWS USERS PLEASE NOTE: Be sure to add the ITK bin directory to your path. This will resolve the ITK dll's at run time.

Building All of the Examples

Many of the examples in the ITK Wiki Examples Collection require VTK. You can build all of the the examples by following these instructions. If you are a new VTK user, you may want to try the Superbuild which will build a proper ITK and VTK.

ItkVtkGlue

ITK >= 4

For examples that use QuickView (which depends on VTK), you must have built ITK with Module_ITKVtkGlue=ON.

ITK < 4

Some of the ITK Examples require VTK to display the images. If you download the entire ITK Wiki Examples Collection, the ItkVtkGlue directory will be included and configured. If you wish to just build a few examples, then you will need to download ItkVtkGlue and build it. When you run cmake it will ask you to specify the location of the ItkVtkGlue binary directory.