CODEWITHSUNDEEP
c++opencvlevenberg-marquardt
Frank-Rene Schäfer
Frank-Rene Schäfer

Reputation: 3352

OpenCV::LMSolver getting a simple example to run

PROBLEM: The documentation for cv::LMSolver at opencv.org is very thin, to say the least. Finding some useful examples on the internet, also, was not possible.

APPROACH: So, I did write some simple code:

#include <opencv2/calib3d.hpp>
#include <iostream>

using namespace cv;
using namespace std;

struct Easy : public LMSolver::Callback {
   Easy() = default;

   virtual bool compute(InputArray f_param, OutputArray f_error, OutputArray f_jacobian) const override
   { 
        Mat  param = f_param.getMat();

        if( f_error.empty() ) f_error.create(1, 1, CV_64F);             // dim(error) = 1
        Mat  error = f_error.getMat();                          

        vector<double> x{param.at<double>(0,0), param.at<double>(1,0)}; // dim(param) = 2
        double         error0 = calc(x);
        error.at<double>(0,0) = error0;

        if( ! f_jacobian.needed() ) return true; 
        else if( f_jacobian.empty() ) f_jacobian.create(1, 2, CV_64F); 
        Mat  jacobian = f_jacobian.getMat();

        double e = 1e-10; // estimate derivatives in epsilon environment       
        jacobian.at<double>(0, 0) = (calc({x[0] + e, x[1]    }) - error0) / e; // d/dx0 (error)
        jacobian.at<double>(0, 1) = (calc({x[0],     x[1] + e}) - error0) / e; // d/dx1 (error)
        return true;
   }

   double calc(const vector<double> x) const { return x[0]*x[0] + x[1]*x[1]; }
};

int main(int argc, char** argv) 
{
    Ptr<Easy>      callback   = makePtr<Easy>();
    Ptr<LMSolver>  solver     = LMSolver::create(callback, 100000, 1e-37);
    Mat            parameters = (Mat_<double>(2,1) << 5, 100);
    solver->run(parameters);
    cout << parameters << endl;
}

QUESTIONS:

Upvotes: 0

Views: 1244

Answers (2)

swk
swk

Reputation: 1

I suppose OP wants to minimize x^2 + y^2 with respect to [x, y].

Because Levenberg-Marquardt method solves a least square problem, the error should be defined as [x, y] so that it minimizes || [x, y] ||^2 = x^2 + y^2.

Another suggestion is that the Jacobian matrix should be provided analytically whenever possible, although this is not crucial in this particular case.

struct Easy : public LMSolver::Callback {
   Easy() = default;

   virtual bool compute(InputArray f_param, OutputArray f_error, OutputArray f_j
acobian) const override
   {
        Mat  param = f_param.getMat();

        if( f_error.empty() ) f_error.create(2, 1, CV_64F);             // dim(error) = 2
        Mat  error = f_error.getMat();

        vector<double> x{param.at<double>(0,0), param.at<double>(1,0)}; // dim(param) = 2
        error.at<double>(0, 0) = x[0];
        error.at<double>(1, 0) = x[1];

        if( ! f_jacobian.needed() ) return true;
        else if( f_jacobian.empty() ) f_jacobian.create(2, 2, CV_64F);
        Mat  jacobian = f_jacobian.getMat();

        jacobian.at<double>(0, 0) = 1;
        jacobian.at<double>(0, 1) = 0;
        jacobian.at<double>(1, 0) = 0;
        jacobian.at<double>(1, 1) = 1;

        return true;
   }
};

Upvotes: 0

user4442671
user4442671

Reputation:

What does the return value of LMSolver::Callback::compute() report to the caller?

Thankfully, opencv is opensource, so we might be able to figure this out simply by checking out the code.

Looking at the source code on Github, I found that all of the calls to compute() look like:

if( !cb->compute(x, r, J) )
    return -1;

Returning false simply causes the solver to bail out. So it seems that the return value of the callback's compute() is simply whether the generation of the jacobian was successful or not.

Currently, it finds the minimum at (-9e-07,4e-5). How can the precision be improved?

If anything, you should at least compare the return value of run() against your maximum iteration count to make sure that it did, in fact, converge as much as it could.

Upvotes: 2

Related Questions