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

[<bl...@us...>]

Revision: 3868
          http://hugin.svn.sourceforge.net/hugin/?rev=3868&view=rev
Author:   blimbo
Date:     2009-05-21 17:07:15 +0000 (Thu, 21 May 2009)

Log Message:
-----------
Test program test_line_dist.cpp to work out curve-point distance, can't quite get it to work yet

Modified Paths:
--------------
    hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/ParameterEstimator.cpp
    hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/ProcessImage.h

Added Paths:
-----------
    hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/test_line_dist.cpp

Modified: hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/ParameterEstimator.cpp
===================================================================
--- hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/ParameterEstimator.cpp	2009-05-20 18:27:37 UTC (rev 3867)
+++ hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/ParameterEstimator.cpp	2009-05-21 17:07:15 UTC (rev 3868)
@@ -1,5 +1,5 @@
-#include <cmath>
-#include "ParameterEstimator.h"
+#include <cmath>
+#include "ParameterEstimator.h"
 #include <iostream>
 #include "vigra/diff2d.hxx"
 #include <boost/numeric/ublas/matrix.hpp>
@@ -10,6 +10,8 @@
 using namespace boost::numeric::ublas;
 using namespace std;
 
+double a,b,c,p,q;
+
 ParameterEstimator::ParameterEstimator(unsigned int m, double delta){
 
 	// Number of data objects required for an exact 
@@ -135,20 +137,21 @@
  */
 
 
-/* model to be fitted to measurements: x_i = p[0]*exp(-p[1]*i) + p[2], i=0...n-1 */
+
+
 void func(double *p, double *x, int m, int n, void *data){
 	register int i;
 
 	std::vector<const vigra::Point2D *> * dat = (std::vector<const vigra::Point2D *> *) data;
 
   	for(i=0; i < n; ++i){
-		//cout << (*dat)[i]->x << "," << (*dat)[i]->y << endl;
-    		//[i]= p[0] * exp(-p[1]*i) + p[2];
+		//std::cout << "x = " << (*dat)[i]->x << ", y =" << (*dat)[i]->y << endl;
+    		//[i]= ???
   	}
 }
 
 void jacfunc(double *p, double *jac, int m, int n, void *data){}
-
+
 void ParameterEstimator::leastSquaresEstimate(std::vector<const vigra::Point2D *> &data, std::vector<double> &parameters)
 {
 
@@ -159,8 +162,9 @@
 	register int i;
 	int ret;
 
-  	//for(i=0; i<n; ++i)
-    	//	x[i]= data[i]->x;
+  	for(i=0; i < n; ++i){
+    		x[i]= 1;
+	}
 
   	/* initial parameters estimate: (1.0, 1.0, 1.0) */
   	p[0]=1.0; p[1]=1.0; p[2]=1.0;
@@ -175,15 +179,17 @@
   	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> &parameters, const vigra::Point2D &data, double x_min, double x_max){
 
+
 	double shortest = 100000000;
-	for (int x = x_min; x < x_max; x++){
+
+	for (int x = (int)x_min; x < x_max; x++){
 		double y = (parameters[0]*(x*x)) + (parameters[1]*x) + parameters[2];	
 		double distancesquared = abs((data->x - x) + (data->y - y));
 		if(distancesquared < shortest){
@@ -191,8 +197,35 @@
 		}
 		//std::cout << x << "\t=\t" << distancesquared << std::endl;
 	}
+	
+	
+	//std::cout << "Brute force:\t" << shortest << std::endl;
+
+
+	if (1){
+
+		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;
+
+
+		shortest = -(b/(2*a))
+				- (((6*a2) - (3*a2*b2) + (12*a3*c) - (12*a3*q))
+				/(3 * pow(2,(double)(2/3)) * 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)))),(double)(1/3))))
+				+ (1/(6*pow(2,(double)(1/3))*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))),(double)(1/3)));
+
+		std::cout << "Algebreic:\t" << shortest << std::endl;
+	}
+
 	return (shortest < this->deltasquared);
-	
-}
+	
+}
 
 

Modified: hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/ProcessImage.h
===================================================================
--- hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/ProcessImage.h	2009-05-20 18:27:37 UTC (rev 3867)
+++ hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/ProcessImage.h	2009-05-21 17:07:15 UTC (rev 3868)
@@ -4,7 +4,7 @@
 #include "vigra/impex.hxx"
 #include "vigra/diff2d.hxx"
 
-void ransac(std::vector<vigra::Point2D>&);
+void ransac(std::vector<vigra::Point2D>&,unsigned int,unsigned int);
 void extract_coords(vigra::BImage&);
 void sub_image(vigra::BImage& image);
 void detect_edge(vigra::BImage&, std::string&, std::string&, std::string&);

Added: hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/test_line_dist.cpp
===================================================================
--- hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/test_line_dist.cpp	                        (rev 0)
+++ hugin/branches/gsoc2009_lenscalibration/src/lens_calibrate/test_line_dist.cpp	2009-05-21 17:07:15 UTC (rev 3868)
@@ -0,0 +1,45 @@
+#include <cmath>
+#include <iostream>
+
+using namespace std;
+
+/*
+
+http://answers.yahoo.com/question/index?qid=20070109172252AAP34wx
+Min distance from (3,2) to y = x^2 + 2 should = 2.23607
+
+*/ 
+
+double p = 3;
+double q = 2;
+double a = 1;
+double b = 0;
+double c = 2;
+double a2 = a*a;
+double a3 = a2*a;
+double a4 = a2*a2;
+double b2 = b*b;
+
+int main(){
+
+	double shortest_squared = 	
+					-(b/(2*a))
+					- (6*a2 - 3*a2*b2 + 12*a3*c - 12*a3*q)
+					/(3*pow(2,(double)(2/3))*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))
+					),(double)(1/3)))
+					+ (1/(6*pow(2,(double)(1/3))*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)))
+					,(double)(1/3));
+
+	std::cout << "Dist^2:\t" << shortest_squared << std::endl;
+	std::cout << "Dist:\t" << sqrt(shortest_squared) << std::endl;
+
+	// closest point is (1,3)
+	std::cout << "Correct:\t" << sqrt ((3-1)*(3-1) + (2 - 3)*( 2 - 3)) << std::endl;
+
+	return (1);
+	
+}


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.