| |
 |
|
|
Science Forum Index » Nonlinear Science Forum » help needed with levmar minimizer
Page 1 of 1
|
| Author |
Message |
| Guest |
 Posted: Wed Aug 23, 2006 7:30 am |
|
|
|
|
Hi,
I am trying to write a program to fit the sphere to measured data,
using the least squares fit. For that purpose I am using a levmar
function (http://www.ics.forth.gr/~lourakis/levmar). The problem is the
documentation, or rather the lack of it. I don't know how to pass my
data to levmar solver, to minimize my objective function. Examples
provided with program deal with analytical problems only, not the real
ones involving data.
I spent past few days trying to figure out how to pass data to levmar
minimizer, so I am pretty desperate. I hope you can help me.. Of course
I can send a code I wrote, so you can see where I got it wrong.
regards
Reply »
p.s.
here is the source code
#include <iostream>
#include <cmath>
#include "levmar-2.1.2/lm.h"
struct mydata
{
double X[864];
};
// Objective function
void func(double *a, double *X, int m, int n, void *data)
{
double *xt (new double[n]);
double *yt (new double[n]);
double *zt (new double[n]);
double *r (new double[n]);
struct mydata *d = (struct mydata *)data;
int j=0;
for (int i=0; i<n; i+=3)
{
xt[j] = d->X[i] - a[0];
yt[j] = d->X[i+1] - a[1];
zt[j] = d->X[i+2] - a[2];
r[j] = std::sqrt(xt[j]*xt[j] + yt[j]*yt[j] + zt[j]*zt[j]);
X[j] = X[j-1]+(r[j] - a[3]);
++j;
}
}
// Jacobian
void jacobian(double *a, double *J, int m, int n, void *data)
{
double *xt (new double[n]);
double *yt (new double[n]);
double *zt (new double[n]);
double *r (new double[n]);
struct mydata *d;
d =(struct mydata *)data;
int j=3;
int k=0;
for (int i=0; i<n; i+=3)
{
xt[k] = d->X[i] - a[0];
yt[k] = d->X[i+1] - a[1];
zt[k] = d->X[i+2] - a[2];
r[k] = std::sqrt(xt[k]*xt[k] + yt[k]*yt[k] + zt[k]*zt[k]);
// Jacobian
J[i] = -xt[k]/r[k];
J[i+1] = -yt[k]/r[k];
J[i+2] = -zt[k]/r[k];
J[j] = -1;
j+=4; ++k;
}
}
int main()
{
int size = 864;
int radius = 5;
const double pi = 3.14159265;
double *X ( new double[size] );
mydata data;
// point data generation, generates two spheres, our fitted sphere
should have a center at (0,0,0) and radius of 5.5
int k = -1;
for(int i=0; i<24; i++)
{
for (int j=0; j<12; j++)
{
if (j%2==0)
{
data.X[++k] = radius*std::sin(j*pi/12)*std::cos(i*pi/12);
data.X[++k] = radius*std::sin(j*pi/12)*std::sin(i*pi/12);
data.X[++k] = radius*std::cos(j*pi/12);
}
else
{
data.X[++k] = (radius+1)*std::sin(j*pi/12)*std::cos(i*pi/12);
data.X[++k] = (radius+1)*std::sin(j*pi/12)*std::sin(i*pi/12);
data.X[++k] = (radius+1)*std::cos(j*pi/12);
}
}
}
for (i=0; i<size; i++)
X[i] = 0;
// Initial estimate
double a[10];
// Position of center point
a[0] = a[1] = a[2] = 0.2;
// radius
a[3] = 4.2;
double *J ( new double[size*4] );
double opts[LM_OPTS_SZ], info[LM_INFO_SZ];
opts[0]=LM_INIT_MU; opts[1]=1E-15; opts[2]=1E-15; opts[3]=1E-20;
opts[4]=LM_DIFF_DELTA;
int ret = dlevmar_der(func, jacobian, a, X, 4, size, 10000, opts,
info, NULL, NULL, (void *)&data);
// Taken from lmdemo.c file
printf("Levenberg-Marquardt returned %d in %g iter, reason
%g\nSolution: ", ret, info[5], info[6]);
for(i=0; i<4; ++i)
printf("%.7g ", a[i]);
printf("\n\nMinimization info:\n");
for(i=0; i<LM_INFO_SZ; ++i)
printf("%g ", info[i]);
printf("\n");
return(0);
} |
|
|
| Back to top |
|
| |
|
Page 1 of 1
All times are GMT - 5 Hours
The time now is Wed Jan 07, 2009 10:39 pm
|
|