From: <bl...@us...> - 2009-05-29 16:36:58
|
Revision: 3892 http://hugin.svn.sourceforge.net/hugin/?rev=3892&view=rev Author: blimbo Date: 2009-05-29 16:36:52 +0000 (Fri, 29 May 2009) Log Message: ----------- Added spline parameter estimator Modified Paths: -------------- hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/SplineEstimator.cpp hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/SplineEstimator.h Modified: hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/SplineEstimator.cpp =================================================================== --- hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/SplineEstimator.cpp 2009-05-29 05:23:39 UTC (rev 3891) +++ hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/SplineEstimator.cpp 2009-05-29 16:36:52 UTC (rev 3892) @@ -1,165 +1,187 @@ #include "SplineEstimator.h" #include <cmath> +#include <iomanip> #include <iostream> #include "vigra/diff2d.hxx" -#include <boost/numeric/ublas/matrix.hpp> -#include <boost/numeric/ublas/io.hpp> #include <levmar/lm.h> #include "MatrixDeterminant.h" #include "Spline.h" -using namespace boost::numeric::ublas; using namespace std; + +double frunge (double x); +double fprunge (double x); +double fpprunge (double x); + double a,b,c,p,q; +bool compare (const vigra::Point2D * i,const vigra::Point2D * j) { + return (i->x<j->x); +} + + ParameterEstimator::ParameterEstimator(unsigned int m, double delta){ - // Number of data objects required for an exact - // estimate (e.g. 2 for a line where the data objects are points) minForEstimate = m; deltasquared = delta*delta; } bool ParameterEstimator::estimate(std::vector<const vigra::Point2D *> &data, std::vector<double> ¶meters){ - parameters.clear(); - if(data.size()<this->minForEstimate) + parameters.clear(); + if(data.size()<this->minForEstimate){ return(false); + } + + sort(data.begin(), data.end(), compare); + if((data[0]->x == data[1]->x)||(data[1]->x == data[2]->x)){ + return(false); + } - /* Quadratic: ax^2 + bx + c = y - Need to estimate a, b and c - - Using Cramer's Rule - http://en.wikipedia.org/wiki/Cramer%27s_rule - [x_1^2, x_1, 1] [a] [y_1] - [x_2^2, x_2, 1]*[b]=[y_2] - [x_3^2, x_3, 1] [c] [y_3] - - Divisions below are using the determinants of - the matrices - - [ y_1 , x_1, 1] - [ y_2 , x_2, 1] - [ y_3 , x_3, 1] - a = --------------- - [x_1^2, x_1, 1] - [x_2^2, x_2, 1] - [x_3^2, x_3, 1] - - [x_1^2, y_1, 1] - [x_2^2, y_2, 1] - [x_3^2, y_3, 1] - b = --------------- - [x_1^2, x_1, 1] - [x_2^2, x_2, 1] - [x_3^2, x_3, 1] - - [x_1^2, x_1, y_1] - [x_2^2, x_2, y_2] - [x_3^2, x_3, y_3] - c = ----------------- - [x_1^2, x_1, 1] - [x_2^2, x_2, 1] - [x_3^2, x_3, 1] + /* + Cubic splines have 2n+2 parameters, where n is the number of + data points. The first n parameters are the x-values. The next + n parameters are the y-values. The last two parameters are + the values of the derivative at the first and last point. For natural + splines, these parameters are unused. + + http://www.extremeoptimization.com/QuickStart/CubicSplinesVB.aspx + */ - matrix<double> m1 (3, 3), m2 (3, 3); - m1 (0,0) = data[0]->x*data[0]->x; - m1 (0,1) = data[0]->x; - m1 (0,2) = 1; - m1 (1,0) = data[1]->x*data[1]->x; - m1 (1,1) = data[1]->x; - m1 (1,2) = 1; - m1 (2,0) = data[2]->x*data[2]->x; - m1 (2,1) = data[2]->x; - m1 (2,2) = 1; + for(int i = 0; i < data.size(); i++){ + parameters.push_back(data[i]->x); + } + for(int i = 0; i < data.size(); i++){ + parameters.push_back(data[i]->y); + } - double d_m1 = determinant(m1); - //std::cout << "Matrix m1:\t" << m1 << std::endl; - //std::cout << "m1 determinant:\t" << d_m1 << std::endl; +// - m2 = m1; - m2 (0,0) = data[0]->y; - m2 (1,0) = data[1]->y; - m2 (2,0) = data[2]->y; + int N = data.size(); + int i; + int ibcbeg; + int ibcend; + int j; + int jhi; + int k; + double t[N]; + double tval; + double y[N]; + double ybcbeg; + double ybcend; + double *ypp; + double yppval; + double ypval; + double yval; - /* - std::cout << "Matrix m2:\t" << m2 << std::endl; - double d_m2_a = determinant(m2); - std::cout << "(a) m2 determinant:\t" << d_m2_a << std::endl; - double a = d_m2_a/d_m1; - */ - double a = determinant(m2)/d_m1; + cout << "SPLINE_CUBIC_SET sets up a cubic spline;\n"; + cout << "SPLINE_CUBIC_VAL evaluates it.\n"; + cout << "\n"; + cout << "Runge's function, evenly spaced knots.\n\n"; + // Points (knots?) are not evenly spaced - is this important? + cout << "I\tT\tY\n\n"; + for ( i = 0; i < N; i++ ) { + t[i] = data[i]->x; + y[i] = data[i]->y; + cout << i << "\t" << t[i] << "\t" << y[i] << "\n"; + } - m2 = m1; - m2 (0,1) = data[0]->y; - m2 (1,1) = data[1]->y; - m2 (2,1) = data[2]->y; - /* - std::cout << "Matrix m2:\t" << m2 << std::endl; - double d_m2_b = determinant(m2); - std::cout << "(b) m2 determinant:\t" << d_m2_b << std::endl; - double b = d_m2_b/d_m1; - */ - double b = determinant(m2)/d_m1; + ibcbeg = 1; + ybcbeg = fprunge ( t[0] ); - m2 = m1; - m2 (0,2) = data[0]->y; - m2 (1,2) = data[1]->y; - m2 (2,2) = data[2]->y; + ibcend = 1; + ybcend = fprunge ( t[N-1] ); - /* - std::cout << "Matrix m2:\t" << m2 << std::endl; - double d_m2_c = determinant(m2); - std::cout << "(c) m2 determinant:\t" << d_m2_c << std::endl; - double c = d_m2_c/d_m1; - */ - double c = determinant(m2)/d_m1; + cout << "\nBoundary condition 1 at both ends:\n"; + cout << "Y'(left) = " << ybcbeg << "\n"; + cout << "Y'(right) = " << ybcend << "\n"; - /* - std::cout << "Parameter estimates:" << std::endl; - std::cout << "a:\t" << a << std::endl; - std::cout << "b:\t" << b << std::endl; - std::cout << "c:\t" << c << std::endl << std::endl; - */ + // Are these values the derivatives needed to complete the paramaters + // or do i get them by calling spline_cubic_val for 1st and last values? + //parameters.push_back(ybcbeg); + //parameters.push_back(ybcend); + // Can we return here? - parameters.push_back(a); - parameters.push_back(b); - parameters.push_back(c); - return(true); + ypp = spline_cubic_set ( N, t, y, ibcbeg, ybcbeg, ibcend, ybcend ); + cout << "\n"; + cout << "SPLINE''(T)\tF''(T):\n"; + cout << "\n"; + for ( i = 0; i < N; i++ ){ + cout << ypp[i] << "\t" << fpprunge (t[i]) << "\n"; + } -} -/*****************************************************************************/ -/* - * Compute the line parameters [n_x,n_y,a_x,a_y] - */ + cout << "\n"; + cout << "T\tSPLINE(T)\tF(T)\n"; + cout << "\n"; + for ( i = 0; i <= N; i++ ){ + + if ( i == 0 ){ + jhi = 1; + }else if ( i < N ){ + jhi = 2; + }else{ + jhi = 2; + } + for ( j = 1; j <= jhi; j++ ){ + if ( i == 0 ){ + tval = t[0] - 1.0; + }else if ( i < N ){ + tval = ( + ( double ) ( jhi - j + 1 ) * t[i-1] + + ( double ) ( j - 1 ) * t[i] ) + / ( double ) ( jhi ); + }else{ + if ( j == 1 ){ + tval = t[N-1]; + }else{ + tval = t[N-1] + 1.0; + } + } + yval = spline_cubic_val ( N, t, tval, y, ypp, &ypval, &yppval ); -void func(double *p, double *x, int m, int n, void *data){ - register int i; + cout << tval << "\t" << yval << "\t\t" << frunge (tval) << "\n"; - std::vector<const vigra::Point2D *> * dat = (std::vector<const vigra::Point2D *> *) data; + // I think these are the derivitive parameters we need + // If so, we don't really need to be in either of these for loops, just call + // spline_cubic_val on 1st and last point + if ((tval == data[0]->x) || (tval == data[N-1]->x)){ + //cout << i << " " << tval << " deriv:" << ypval << endl; + parameters.push_back(ypval); + } + + } + + + } + + delete [] ypp; - for(i=0; i < n; ++i){ - //std::cout << "x = " << (*dat)[i]->x << ", y =" << (*dat)[i]->y << endl; - //[i]= ??? + cout << endl << "Spline Parameters:" << endl << endl; + for ( i = 0; i < parameters.size(); i++ ) { + cout << i << ":\t" << parameters[i] << endl; } + + cout << endl << "====================================" << endl << endl; + + return(true); + } -void jacfunc(double *p, double *jac, int m, int n, void *data){} +void ParameterEstimator::leastSquaresEstimate(std::vector<const vigra::Point2D *> &data, std::vector<double> ¶meters){ -void ParameterEstimator::leastSquaresEstimate(std::vector<const vigra::Point2D *> &data, std::vector<double> ¶meters) -{ + std::cout << "Least squares estimate..." << std::endl; - std::cout << "Least squares estimate using LEVMAR..." << std::endl; - const int n= data.size(), m = 3; // n measurements, 3 parameters double p[m], x[n], opts[LM_OPTS_SZ], info[LM_INFO_SZ]; + + /* register int i; int ret; @@ -167,18 +189,19 @@ x[i]= 1; } - /* initial parameters estimate: (1.0, 1.0, 1.0) */ + // initial parameters estimate: (1.0, 1.0, 1.0) p[0]=1.0; p[1]=1.0; p[2]=1.0; - /* optimization control parameters; passing to levmar NULL instead of opts reverts to defaults */ + // optimization control parameters; passing to levmar NULL instead of opts reverts to defaults opts[0]=LM_INIT_MU; opts[1]=1E-15; opts[2]=1E-15; opts[3]=1E-20; opts[4]=LM_DIFF_DELTA; // relevant only if the finite difference Jacobian version is used - /* invoke the optimization function */ + // invoke the optimization function ret=dlevmar_der(func, jacfunc, p, x, m, n, 1000, opts, info, NULL, NULL, &data); // with analytic Jacobian printf("Levenberg-Marquardt returned in %g iter, reason %g, sumsq %g [%g]\n", info[5], info[6], info[1], info[0]); printf("Best fit parameters: %.7g %.7g %.7g\n", p[0], p[1], p[2]); + */ parameters.push_back(p[0]); parameters.push_back(p[1]); @@ -217,4 +240,94 @@ } +//****************************************************************************80 +double frunge (double x){ + +//****************************************************************************80 +// +// Purpose: +// +// FRUNGE sets the Runge data values. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 13 January 2007 +// +// Author: +// +// John Burkardt +// + + double fx; + + fx = 1.0 / ( 1.0 + 25.0 * x * x ); + + return fx; +} +//****************************************************************************80 + +double fprunge(double x){ + +//****************************************************************************80 +// +// Purpose: +// +// FPRUNGE sets the Runge derivative values at the endpoints. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 13 January 2007 +// +// Author: +// +// John Burkardt +// + + double bot; + double fx; + + bot = 1.0 + 25.0 * x * x; + fx = -50.0 * x / ( bot * bot ); + + return fx; +} +//****************************************************************************80 + +double fpprunge (double x){ + +//****************************************************************************80 +// +// Purpose: +// +// FPPRUNGE sets the Runge second derivative values at the endpoints. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 13 January 2007 +// +// Author: +// +// John Burkardt +// + + double bot; + double fx; + + bot = 1.0 + 25.0 * x * x; + fx = ( -50.0 + 3750.0 * x * x ) / ( bot * bot * bot ); + + return fx; +} Modified: hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/SplineEstimator.h =================================================================== --- hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/SplineEstimator.h 2009-05-29 05:23:39 UTC (rev 3891) +++ hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/SplineEstimator.h 2009-05-29 16:36:52 UTC (rev 3892) @@ -58,7 +58,7 @@ protected: unsigned int minForEstimate; -private: +private: double deltasquared; //given line L and point P, if dist(L,P)^2 < m_delta^2 then the point is on the line }; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |