From: <bl...@us...> - 2009-06-02 00:04:00
|
Revision: 3903 http://hugin.svn.sourceforge.net/hugin/?rev=3903&view=rev Author: blimbo Date: 2009-06-02 00:03:55 +0000 (Tue, 02 Jun 2009) Log Message: ----------- Fixed array bug, agree function now measures distance from point to straght line between 2 spline points Modified Paths: -------------- hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/ProcessImage.cpp hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/Spline.cpp 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/ProcessImage.cpp =================================================================== --- hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/ProcessImage.cpp 2009-06-01 20:38:48 UTC (rev 3902) +++ hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/ProcessImage.cpp 2009-06-02 00:03:55 UTC (rev 3903) @@ -36,15 +36,18 @@ cout << coords.size() << " 'edge' coordinates found." << endl; vector<double> curveParameters; - double desiredProbabilityForNoOutliers = 0.9; + double desiredProbabilityForNoOutliers = 0.9; double maximalOutlierPercentage = 0.8; - // Use 3 points to calculate a,b and c parameters - // Point must be within 0.2 units of curve to agree - ParameterEstimator curveEstimator(3, 1); + // Use 4 points + // Point must be within 10 units of curve to agree + ParameterEstimator curveEstimator(4, 0.5); cout << endl << "Computing RANSAC..." << endl; double usedData = Ransac::compute(curveParameters, &curveEstimator, coords, desiredProbabilityForNoOutliers, maximalOutlierPercentage); - cout << "y = " << curveParameters[0] << "x^2 + " << curveParameters[1] << "x + " << curveParameters[2] << endl; + for (int i = 0; i < curveParameters.size(); i++){ + cout << "Parameter " << i << ":\t" << curveParameters[i] << endl; + } + cout << "Percentage of points which were used for final estimate: "<< usedData/coords.size() << " % (" << usedData << ")" << endl << endl; exit(1); Modified: hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/Spline.cpp =================================================================== --- hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/Spline.cpp 2009-06-01 20:38:48 UTC (rev 3902) +++ hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/Spline.cpp 2009-06-02 00:03:55 UTC (rev 3903) @@ -9,6 +9,7 @@ # include <iomanip> # include <cmath> # include <ctime> +#include <string.h> using namespace std; Modified: hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/SplineEstimator.cpp =================================================================== --- hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/SplineEstimator.cpp 2009-06-01 20:38:48 UTC (rev 3902) +++ hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/SplineEstimator.cpp 2009-06-02 00:03:55 UTC (rev 3903) @@ -9,22 +9,24 @@ using namespace std; - double frunge (double x); double fprunge (double x); double fpprunge (double x); +int dot(const vigra::Point2D A, const vigra::Point2D B, const vigra::Point2D C); +int cross(const vigra::Point2D A, const vigra::Point2D B, const vigra::Point2D C); +double distance(const vigra::Point2D A, const vigra::Point2D B); +double linePointDist(const vigra::Point2D A, const vigra::Point2D B, const vigra::Point2D C, bool isSegment); -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){ minForEstimate = m; - deltasquared = delta*delta; + distance = delta; + distancesquared = delta*delta; } bool ParameterEstimator::estimate(std::vector<const vigra::Point2D *> &data, std::vector<double> ¶meters){ @@ -35,8 +37,10 @@ } sort(data.begin(), data.end(), compare); - if((data[0]->x == data[1]->x)||(data[1]->x == data[2]->x)){ - return(false); + for (int t = 0; t < data.size()-1; t++){ + if (data[t]->x == data[t+1]->x){ + return(false); + } } /* @@ -66,9 +70,9 @@ int j; int jhi; int k; - double t[N]; + double tval; - double y[N]; + double ybcbeg; double ybcend; double *ypp; @@ -76,6 +80,12 @@ double ypval; double yval; + //double t[N]; + //double y[N]; + + double * t = new double[N]; + double * y = new double[N]; + cout << "SPLINE_CUBIC_SET sets up a cubic spline;\n"; cout << "SPLINE_CUBIC_VAL evaluates it.\n"; cout << "\n"; @@ -170,6 +180,9 @@ cout << endl << "====================================" << endl << endl; + delete [] t; + delete [] y; + return(true); } @@ -178,65 +191,40 @@ std::cout << "Least squares estimate..." << 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; - - for(i=0; i < n; ++i){ - x[i]= 1; + for (int i = 0; i < data.size(); i++){ + cout << data[i]->x << "," << data[i]->y << endl; } - // 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 - 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 - 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]); - parameters.push_back(p[2]); + //parameters.push_back(p[0]); + //parameters.push_back(p[1]); + //parameters.push_back(p[2]); } bool ParameterEstimator::agree(std::vector<double> ¶meters, const vigra::Point2D &data){ - double p = data->x; - double q = data->y; - double a = parameters[0]; - double b = parameters[1]; - double c = parameters[2]; - double a2 = a*a; - double a3 = a2*a; - double a4 = a2*a2; - double b2 = b*b; + //cout << endl << "Agree Parameters:" << endl; + //for (int i = 0; i < parameters.size()-2; i++ ){ + // cout << i << ":\t" << parameters[i] << endl; + //} - double x = - -(b/(2*a)) - - (6*a2 - 3*a2*b2 + 12*a3*c - 12*a3*q) - /(3*pow(2,(2.0/3.0))*a2*pow((54*a3*b + 108*a4*p + - sqrt(pow((54*a3*b + 108*a4*p),2) + 4*pow((6*a2 - 3*a2*b2 + 12*a3*c - 12*a3*q),3)) - ),(1.0/3.0))) - + (1/(6*pow(2,(1.0/3.0))*a2)) - * pow( - (54*a3*b + 108*a4*p + sqrt(pow((54*a3*b + 108*a4*p),2) - + 4*pow((6*a2 - 3*a2*b2 + 12*a3*c - 12*a3*q),3))) - ,(1.0/3.0)); - - double y = (a*(x*x)) + (b*x) + c; + vector<vigra::Point2D> coords; + double min_distance = 100000; + + for (int i = 0; i < this->minForEstimate; i++ ){ + //cout << "pushing " << i << " and " << i+this->minForEstimate << endl; + coords.push_back(vigra::Point2D(parameters[i],parameters[i+this->minForEstimate])); + } - double distancesquared = abs((p - x) + (q - y)); + for (int i = 0; i < coords.size()-1; i++ ){ + double distance = linePointDist(coords[i],coords[i+1],data,true); + if(distance < min_distance)min_distance = distance; + //cout << "Distance = " << distance << endl; + } - return (distancesquared < this->deltasquared); + //cout << "Min distance = " << min_distance << endl; + + return(min_distance < this->distance); } @@ -331,3 +319,52 @@ return fx; } + + +// Compute the dot product AB . BC +int dot(const vigra::Point2D A, const vigra::Point2D B, const vigra::Point2D C){ + + int AB[2]; + int BC[2]; + AB[0] = B->x-A->x; + AB[1] = B->y-A->y; + BC[0] = C->x-B->x; + BC[1] = C->y-B->y; + int dot = AB[0] * BC[0] + AB[1] * BC[1]; + return dot; + +} +// Compute the cross product AB x AC +int cross(const vigra::Point2D A, const vigra::Point2D B, const vigra::Point2D C){ + + int AB[2]; + int AC[2]; + AB[0] = B->x-A->x; + AB[1] = B->y-A->y; + AC[0] = C->x-A->x; + AC[1] = C->y-A->y; + int cross = AB[0] * AC[1] - AB[1] * AC[0]; + return cross; +} + + +// Compute the distance from A to B +double distance(const vigra::Point2D A, const vigra::Point2D B){ + int d1 = A->x - B->x; + int d2 = A->y - B->y; + return sqrt(d1*d1+d2*d2); +} + +// Compute the distance from AB to C +// If isSegment is true, AB is a segment, not a line. +double linePointDist(const vigra::Point2D A, const vigra::Point2D B, const vigra::Point2D C, bool isSegment){ + + double dist = cross(A,B,C) / distance(A,B); + if(isSegment){ + int dot1 = dot(A,B,C); + if(dot1 > 0)return distance(B,C); + int dot2 = dot(B,A,C); + if(dot2 > 0)return distance(A,C); + } + return abs(dist); +} Modified: hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/SplineEstimator.h =================================================================== --- hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/SplineEstimator.h 2009-06-01 20:38:48 UTC (rev 3902) +++ hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/SplineEstimator.h 2009-06-02 00:03:55 UTC (rev 3903) @@ -59,7 +59,7 @@ unsigned int minForEstimate; private: - double deltasquared; //given line L and point P, if dist(L,P)^2 < m_delta^2 then the point is on the line + double distance,distancesquared; //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. |