## [Hugin-cvs] SF.net SVN: hugin:[3849] hugin/branches/gsoc2009_lenscalibration/src/ lens_calibrate

 [Hugin-cvs] SF.net SVN: hugin:[3849] hugin/branches/gsoc2009_lenscalibration/src/ lens_calibrate From: - 2009-05-12 22:19:49 ```Revision: 3849 http://hugin.svn.sourceforge.net/hugin/?rev=3849&view=rev Author: blimbo Date: 2009-05-12 22:19:39 +0000 (Tue, 12 May 2009) Log Message: ----------- Written RANSAC curve estimator using quadratic model Modified Paths: -------------- hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/ParameterEstimator.cpp Added Paths: ----------- hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/MatrixDeterminant.h Added: hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/MatrixDeterminant.h =================================================================== --- hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/MatrixDeterminant.h (rev 0) +++ hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/MatrixDeterminant.h 2009-05-12 22:19:39 UTC (rev 3849) @@ -0,0 +1,33 @@ +#ifndef MATRIX_DETERMINANT_H +#define MATRIX_DETERMINANT_H + +#include + +namespace ublas = boost::numeric::ublas; + +template +double determinant(ublas::matrix_expression const& mat_r){ + + double det = 1.0; + + matrix_T mLu(mat_r() ); + ublas::permutation_matrix pivots(mat_r().size1() ); + + int is_singular = lu_factorize(mLu, pivots); + + if(!is_singular){ + + for (std::size_t i=0; i < pivots.size(); ++i){ + if (pivots(i) != i){ + det *= -1.0; + } + det *= mLu(i,i); + } + }else{ + det = 0.0; + } + + return det; +} + +#endif Modified: hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/ParameterEstimator.cpp =================================================================== --- hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/ParameterEstimator.cpp 2009-05-11 22:56:37 UTC (rev 3848) +++ hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/ParameterEstimator.cpp 2009-05-12 22:19:39 UTC (rev 3849) @@ -2,7 +2,12 @@ #include "ParameterEstimator.h" #include #include "vigra/diff2d.hxx" +#include +#include +#include "MatrixDeterminant.h" +using namespace boost::numeric::ublas; + ParameterEstimator::ParameterEstimator(unsigned int m){ // Number of data objects required for an exact @@ -12,43 +17,101 @@ bool ParameterEstimator::estimate(std::vector &data, std::vector ¶meters){ - // Quadratic y = ax2 + bx + c + parameters.clear(); + if(data.size()minForEstimate) + 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] + // (m1) + // + // [ 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] - // When x == 0. Ok to take average of 3 c's? - double c = (double)(data[0]->y + data[1]->y + data[2]->y)/3; - // Viète's formulas. Ok to take average of 2 b/c's? - double a = ( (c/(data[0]->x * data[1]->x)) + (c/(data[1]->x * data[2]->x)) )/2; - double b = ( (-1*(data[0]->x + data[1]->x)*a) + (-1*(data[1]->x + data[2]->x)) )/2; + matrix 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 (unsigned i = 0; i < data.size(); i++){ - // std::cout << "Estimator: " << i << " " << data[i]->x << "," << data[i]->y << std::endl; - //} + double d_m1 = determinant(m1); + std::cout << "Matrix m1:\t" << m1 << std::endl; + std::cout << "m1 determinant:\t" << d_m1 << std::endl; - parameters.clear(); - if(data.size()minForEstimate) - return(0); - //double nx = data[1]->y - data[0]->y; - //double ny = data[0]->x - data[1]->x; - //double norm = sqrt(nx*nx + ny*ny); - - //parameters.push_back(nx/norm); - //parameters.push_back(ny/norm); - //parameters.push_back(data[0]->x); - //parameters.push_back(data[0]->y); + m2 = m1; + m2 (0,0) = data[0]->y; + m2 (1,0) = data[0]->y; + m2 (2,0) = data[0]->y; + 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; + + m2 = m1; + m2 (0,1) = data[1]->y; + m2 (1,1) = data[1]->y; + m2 (2,1) = data[1]->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; + + m2 = m1; + m2 (0,2) = data[2]->y; + m2 (1,2) = data[2]->y; + m2 (2,2) = data[2]->y; + + 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; + + 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; + parameters.push_back(a); parameters.push_back(b); parameters.push_back(c); + return(true); - char params[3] = {'a', 'b', 'c'}; - for (unsigned i = 0; i < parameters.size(); i++){ - std::cout << "RANSAC curve parameter: " << params[i] << " " << parameters[i] << std::endl; - } - - return(1); } /*****************************************************************************/ /* This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```