itkExponentialDecayFunction.h

Go to the documentation of this file.

/*=============================================================================

  NifTK: A software platform for medical image computing.

  Copyright (c) University College London (UCL). All rights reserved.

  This software is distributed WITHOUT ANY WARRANTY; without even
  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  PURPOSE.

  See LICENSE.txt in the top level directory for details.

=============================================================================*/

#ifndef itkExponentialDecayFunction_h
#define itkExponentialDecayFunction_h

#include <vnl/vnl_least_squares_function.h>


class ITK_EXPORT ExponentialDecayFunction : public vnl_least_squares_function
{
public:

ExponentialDecayFunction(unsigned int nDataPoints,
       vnl_vector< double > &xData,
       vnl_vector< double > &yData,
       bool flgUseGradient ) 
  : vnl_least_squares_function(2, nDataPoints,
         flgUseGradient ? use_gradient : no_gradient) 
  {
    m_nDataPoints = nDataPoints;

    m_xData = &xData;
    m_yData = &yData;
  }

  double compute(unsigned int i, double k, double s) 
  {
    return compute( (*m_xData)[i], k, s );
  }

  double compute(double x, double k, double s) 
  {
    return k*exp( -s*x );
  }

  double compute_dydk(unsigned int i, double k, double s) 
  {
    return compute_dydk( (*m_xData)[i], k, s );
  }

  double compute_dydk(double x, double k, double s) 
  {
    return exp( -s*x );
  }

  double compute_dydsig(unsigned int i, double k, double s) 
  {
    return compute_dydsig( (*m_xData)[i], k, s );
  }

  double compute_dydsig(double x, double k, double s) 
  {
    return -k*x*exp( -s*x );
  }


  void f(const vnl_vector<double> &a, vnl_vector<double> &residual) override 
  {
    for (unsigned int i=0; i<m_nDataPoints; ++i)       
    {
      residual[i] = compute(i, a[0], a[1] ) - (*m_yData)[i];
    }
  }


  void gradf(const vnl_vector<double> &a, vnl_matrix<double> &J) override 
  {
    for (unsigned int i=0; i<m_nDataPoints; ++i) {
      J(i,0) = compute_dydk(i, a[0], a[1] );
    }

    for (unsigned int i=0; i<m_nDataPoints; ++i) {
      J(i,1) = compute_dydsig(i, a[0], a[1] );
    }
  }

protected:

  unsigned int m_nDataPoints;

  vnl_vector< double > *m_xData;
  vnl_vector< double > *m_yData;
  

};

#endif