[Insight-users] vnl_powell: three parameter fit

[Mathieu Malaterre]

Hello,

  Is there anyone here using vnl_powell (*) ? I am trying to estimate
the 3 parameters of a fit function, on let say 15 estimated signal
value. The function I need to pass to the powell optimizer is:

      || F(a,b,c) - Signal ||²

It this correct ? Are there any other restriction on the function to minimize ?

Thanks,
-Mathieu
(*) I know the vnl version is patented, but the idea should remain the
same with another powell implementation


#include <vnl/vnl_vector.h>
#include <vnl/vnl_double_2.h>
#include <vnl/algo/vnl_powell.h>
#include <vnl/vnl_cost_function.h>

#include <iostream>

class dummy_func : public vnl_cost_function
{
public:

  dummy_func(int n) : vnl_cost_function(n) {}

  static const unsigned int nsignals = 15;
  static double compute(double t, double a, double b, double c)
    {
    return a * ( 1. - b * exp( -t / c ) );
    }
  double f(const vnl_vector<double>& x)
  {
    assert((int)x.size()==dim);

    vnl_vector<double> computed(dim);
    for (unsigned int i=0; i<nsignals; ++i)
      computed[i] = compute(i, x(0), x(1), x(2) );

    static const double signal[] = {
      -13.5942,
      -17.7883,
      -23.0262,
      -29.5675,
      -37.7366,
      -47.9385,
      -60.6792,
      -76.5904,
      -96.4611,
      -121.277,
      -152.267,
      -190.97,
      -239.304,
      -299.666,
      -375.049,
    };
    double sum = 0;
    for (unsigned int i=0; i<nsignals; ++i)
      {
      sum += (computed[i]-signal[i])*(computed[i]-signal[i]);
      }
    return sum;
  }
};

int main(int, char *[])
{
//  for(unsigned int i = 0; i < dummy_func::nsignals; ++i)
//    {
//    std::cout << dummy_func::compute(i, 3.26, 5.17, -4.5) << std::endl;
//    }

  const unsigned n = 3;

  vnl_vector<double> x(n);
  x.fill(1);
  dummy_func cost1(n);
  vnl_powell powell(&cost1);
  powell.minimize(x);

  std::cout << x << std::endl;

  return 0;
}