LevenbergMarquardtBase.cpp

/* 
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include "LevenbergMarquardtBase.hh"

using namespace BIAS;

int LevenbergMarquardtBase::
LM_Compute(const int numResiduals, const int numParameters,
           BIAS::Vector<double> &startX, BIAS::Vector<double> &resultX,
           BIAS::Vector<double> &residuals)
{
  // Call Levenberg-Marquardt algorithm from Minpack wrapper
  int res = 0;
  LM_squaredSum_ = 0;
  LM_numIters_ = 0;
  LM_parameters_.newsize(numParameters);
  LM_residuals_.newsize(numResiduals);
  LM_Jacobian_.newsize(numResiduals, numParameters);
  res = BIAS::LevenbergMarquardtExtended(&g_MinpackFunction, this,
                                          numResiduals, numParameters,
                                          startX, resultX,
                                          LM_tolerance_, LM_maxIters_);

  // Evaluate results
  if (res >= 0) {
    residuals = LM_residuals_;
    LM_squaredSum_ = 0;
    for (int j = 0; j < numResiduals; j++)
      LM_squaredSum_ += LM_residuals_[j]*LM_residuals_[j];
  }
  return res;
}

int LevenbergMarquardtBase::
LM_ComputeWithoutJacobian(const int numResiduals,
                          const int numParameters,
                          BIAS::Vector<double> &startX,
                          BIAS::Vector<double> &resultX,
                          BIAS::Vector<double> &residuals,
                          const double epsilon)
{
  // Call Levenberg-Marquardt algorithm from Minpack wrapper
  int res = 0;
  LM_squaredSum_ = 0;
  LM_numIters_ = 0;
  LM_parameters_.newsize(numParameters);
  LM_residuals_.newsize(numResiduals);
  LM_Jacobian_.newsize(numResiduals, numParameters);
  LM_Jacobian_.SetZero();
  res = BIAS::LevenbergMarquardtExtended(&g_MinpackFunctionWithoutJacobian,
                                          this, numResiduals, numParameters,
                                          startX, resultX, epsilon,
                                          LM_tolerance_, LM_maxIters_);

  // Evaluate results
  if (res >= 0) {
    residuals = LM_residuals_;
    LM_squaredSum_ = 0;
    for (int j = 0; j < numResiduals; j++)
      LM_squaredSum_ += LM_residuals_[j]*LM_residuals_[j];
  }

  // Approximate final Jacobian of f
  if (res >= 0) {
    res = BIAS::ComputeJacobian(&g_MinpackFunctionWithoutJacobian, this,
                                numResiduals, resultX, LM_Jacobian_, epsilon);
  } else {
    LM_Jacobian_.SetZero(); // approximation failed!
  }
  return res;
}