Help with curve fitting pls

Beginner Questions
1

Does anyone know how to fit a curve to a series of points with itk?

I currently have a cost function and a Levenberg-Marquad optimizer, but I have no idea how to get the cost function to work - I can’t find documentation anywhere.

I have a series of values (densities) at specific times, and I am trying to fit a gamma-variate function to these points. The gamma-variate function looks like this:

Screenshot_20210301_194234

i.e. I find the best values for A, alpha and beta.

can anyone point me in the right direction? thank you so much!

2

Based on the description, ITK is not really meant to solve this type of problem. I usually use Eigen for this type of problem. If you share some source code, it might ring a bell for someone.

3

thanks for the reply mate -

this is the cost function I’m using. I’ve copied it from an open source program here (GitHub - QIICR/DSC_Analysis: Analysis of Dynamic Susceptibility Contrast (DSC) MRI).


class LMCostFunction: public itk::MultipleValuedCostFunction
{
public:
  typedef LMCostFunction                    Self;
  typedef itk::MultipleValuedCostFunction   Superclass;
  typedef itk::SmartPointer<Self>           Pointer;
  typedef itk::SmartPointer<const Self>     ConstPointer;
  itkNewMacro( Self );
        
  enum { SpaceDimension =  2 };
  unsigned int RangeDimension; 

  enum ModelType { TOFTS_2_PARAMETER = 1 };
        
  typedef Superclass::ParametersType              ParametersType;
  typedef Superclass::DerivativeType              DerivativeType;
  typedef Superclass::MeasureType                 MeasureType, ArrayType;
  typedef Superclass::ParametersValueType         ValueType;
		      
        
  int m_ModelType;

  LMCostFunction()
  {
  }
        
  void SetModelType (int model) {
    m_ModelType = model;
  }
        
  void SetNumberOfValues(unsigned int NumberOfValues)
  {
    RangeDimension = NumberOfValues;
  }
    
  //Cv = measured values 
  void SetCv (const float* cv, int sz) //Self signal Y
  {    
    Cv.set_size (sz);
    for (int i = 0; i < sz; ++i)
      Cv[i] = cv[i];
    //std::cout << "Cv: " << Cv << std::endl;
  }

  
  //Time
  void SetTime (const float* cx, int sz) 
  {
    Time.set_size (sz);
    for( int i = 0; i < sz; ++i )
      Time[i] = cx[i];
    //std::cout << "Time: " << Time << std::endl;
  }
        
  MeasureType GetValue( const ParametersType & parameters) const
  {
    MeasureType measure(RangeDimension);

    ValueType alpha = parameters[0];
    ValueType beta = parameters[1];
            

    ValueType deltaT = Time(1) - Time(0);
    
    if(m_ModelType == TOFTS_2_PARAMETER)
      {
      measure = Cv -(); // function to 'fit' goes here I think 
      }
            
    return measure; 
  }
  
  MeasureType GetFittedFunction( const ParametersType & parameters) const
  {
    MeasureType measure(RangeDimension);

    ValueType alpha = parameters[0];
    ValueType beta = parameters[1];
            

    ValueType deltaT = Time(1) - Time(0);
    
    if(m_ModelType == TOFTS_2_PARAMETER)
      {
      measure = (); //function to 'fit' goes here I think 
      }
            
    return measure; 
  }
    
  //Not going to be used
  void GetDerivative( const ParametersType & /* parameters*/,
                      DerivativeType  & /*derivative*/ ) const
  {   
  }
        
  unsigned int GetNumberOfParameters(void) const
  {
    if(m_ModelType == TOFTS_2_PARAMETER)
      {
      return 2;
      }
  }
        
  unsigned int GetNumberOfValues(void) const
  {
    return RangeDimension;
  }
        
protected:
  virtual ~LMCostFunction(){}
private:
        
  ArrayType Cv, Time; 
        
  ArrayType Convolution(ArrayType X, ArrayType Y) const
  {
    ArrayType Z;
    Z = vnl_convolve(X,Y).extract(X.size(),0);
    return Z;
  };
        
  ArrayType Exponential(ArrayType X) const
  {
    ArrayType Z;
    Z.set_size(X.size());
    for (unsigned int i=0; i<X.size(); i++)
      {
      Z[i] = exp(X(i));
      }
    return Z;
  };
        
  int constraintFunc(ValueType x) const
  {
    if (x<0||x>1)
      return 1;
    else
      return 0;
            
  };
        
        
};

and I’m trying to use it with a Levenberg-Marquad optimizer like this:

itk::LevenbergMarquardtOptimizer::Pointer optimizer = itk::LevenbergMarquardtOptimizer::New();
LMCostFunction::Pointer   costFunction = LMCostFunction::New();

I don’t know how to set up the cost, or ‘measure’ part correctly (I think GetValue() does this). It would have to be something like measure = Cv - (pow(Time, alpha)*exp(-Time/beta)) but it gives me an error that I can’t figure out how to fix.

There is an example of how to use the optimizer here:

https://itk.org/ITKExamples/src/Numerics/Optimizers/LevenbergMarquardtOptimization/Documentation.html

but no example on what the cost function has to look like.

Thanks for your time!