Why i did this silly change

This commit is contained in:
Ankit Gupta 2007-05-14 08:24:35 +00:00
parent 645aac9c61
commit 92ad068512
5 changed files with 335 additions and 1 deletions

2
.gitattributes vendored
View File

@ -1877,6 +1877,8 @@ Principal_component_analysis/demo/Principal_component_analysis/windows/3d/res/id
Principal_component_analysis/doc_tex/Principal_component_analysis/fit.eps -text svneol=unset#application/postscript Principal_component_analysis/doc_tex/Principal_component_analysis/fit.eps -text svneol=unset#application/postscript
Principal_component_analysis/doc_tex/Principal_component_analysis/fit.png -text svneol=unset#image/png Principal_component_analysis/doc_tex/Principal_component_analysis/fit.png -text svneol=unset#image/png
Principal_component_analysis/doc_tex/Principal_component_analysis/teaserLeastSquaresFitting.png -text Principal_component_analysis/doc_tex/Principal_component_analysis/teaserLeastSquaresFitting.png -text
Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments.h -text
Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles.h -text
QP_solver/doc_tex/QP_solver/PkgDescription.tex -text QP_solver/doc_tex/QP_solver/PkgDescription.tex -text
QP_solver/doc_tex/QP_solver/closest_point.eps -text svneol=unset#application/postscript QP_solver/doc_tex/QP_solver/closest_point.eps -text svneol=unset#application/postscript
QP_solver/doc_tex/QP_solver/closest_point.fig -text svneol=unset#application/octet-stream QP_solver/doc_tex/QP_solver/closest_point.fig -text svneol=unset#application/octet-stream

View File

@ -88,6 +88,39 @@ centroid(InputIterator begin,
return ORIGIN + v / (FT)nb_pts; return ORIGIN + v / (FT)nb_pts;
}// end centroid of a 3D point set }// end centroid of a 3D point set
// computes the centroid of a 2D segment set
// takes an iterator range over K::Segment_2
template < typename InputIterator,
typename K >
typename K::Point_2
centroid(InputIterator begin,
InputIterator end,
const K& ,
const typename K::Segment_2*)
{
typedef typename K::FT FT;
typedef typename K::Vector_2 Vector;
typedef typename K::Point_2 Point;
typedef typename K::Segment_2 Segment;
CGAL_precondition(begin != end);
Vector v = NULL_VECTOR;
FT sum_lengths = 0;
for(InputIterator it = begin;
it != end;
it++)
{
const Segment& s = *it;
FT length = std::sqrt(std::abs(s.squared_length()));
// Point c = K().construct_centroid_2_object()(s[0],s[1]);??
v = v + length * (((s[0] - ORIGIN) + (s[1] - ORIGIN))/2);
sum_lengths += length;
}
CGAL_assertion(sum_lengths != 0.0);
return ORIGIN + v / sum_lengths;
} // end centroid of a 2D triangle set
// computes the centroid of a 2D triangle set // computes the centroid of a 2D triangle set
// takes an iterator range over K::Triangle_2 // takes an iterator range over K::Triangle_2
template < typename InputIterator, template < typename InputIterator,
@ -122,7 +155,7 @@ centroid(InputIterator begin,
} // end centroid of a 2D triangle set } // end centroid of a 2D triangle set
// computes the centroid of a 3D triangle set // computes the centroid of a 3D triangle set
// takes an iterator range over K::Triangle_2 // takes an iterator range over K::Triangle_3
template < typename InputIterator, template < typename InputIterator,
typename K > typename K >
typename K::Point_3 typename K::Point_3

View File

@ -23,10 +23,17 @@
#include <CGAL/Object.h> #include <CGAL/Object.h>
#include <CGAL/centroid.h> #include <CGAL/centroid.h>
#include <CGAL/eigen_2.h> #include <CGAL/eigen_2.h>
#include <CGAL/eigen.h>
#include <CGAL/Linear_algebraCd.h>
//#include <CGAL/Kernel_d/Matrix__.h>
#include <CGAL/linear_least_squares_fitting_triangles.h>
#include <CGAL/linear_least_squares_fitting_segments.h>
#include <iterator> #include <iterator>
#include <list> #include <list>
#include <string> #include <string>
#include <vector>
#include <cmath>
CGAL_BEGIN_NAMESPACE CGAL_BEGIN_NAMESPACE

View File

@ -0,0 +1,138 @@
// Copyright (c) 2005 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you may redistribute it under
// the terms of the Q Public License version 1.0.
// See the file LICENSE.QPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL: svn+ssh://gankit@scm.gforge.inria.fr/svn/cgal/trunk/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles.h $
// $Id: linear_least_squares_fitting_2.h 37882 2007-04-03 15:15:30Z spion $
//
// Author(s) : Pierre Alliez and Sylvain Pion and Ankit Gupta
#ifndef CGAL_LINEAR_LEAST_SQUARES_FITTING_SEGMENTS_H
#define CGAL_LINEAR_LEAST_SQUARES_FITTING_SEGMENTS_H
#include <CGAL/basic.h>
#include <CGAL/Object.h>
#include <CGAL/centroid.h>
#include <CGAL/eigen_2.h>
#include <CGAL/eigen.h>
#include <CGAL/Linear_algebraCd.h>
#include <iterator>
#include <vector>
#include <cmath>
CGAL_BEGIN_NAMESPACE
namespace CGALi {
// Fits a line to a 2D segment set.
// Returns a fitting quality (1 - lambda_min/lambda_max):
// 1 is best (zero variance orthogonally to the fitting line);
// 0 is worst (isotropic case, returns a line with horizontal
// direction by default)
template < typename InputIterator, typename K >
typename K::FT
linear_least_squares_fitting_2(InputIterator first,
InputIterator beyond,
typename K::Line_2& line, // best fit line
typename K::Point_2& c, // centroid
const K&, // kernel
const typename K::Segment_2*) // used for indirection
{
// types
typedef typename K::FT FT;
typedef typename K::Line_2 Line;
typedef typename K::Point_2 Point;
typedef typename K::Vector_2 Vector;
typedef typename K::Segment_2 Segment;
typedef typename CGAL::Linear_algebraCd<FT> LA;
typedef typename LA::Matrix Matrix;
// precondition: at least one element in the container.
CGAL_precondition(first != beyond);
// compute centroid
c = centroid(first,beyond,K());
// assemble covariance matrix as a semi-definite matrix.
// Matrix numbering:
// 0
// 1 2
//Final combined covariance matrix for all triangles and their combined mass
FT mass=0.0,cov[3]={0.0,0.0,0.0};
// step 1: assemble the 2nd order moment about the origin.
FT cov_temp[4] = {1.0,0.5,0.5,1.0};
Matrix covariance = (std::sqrt(2)/3.0) * init_Matrix<K>(2,cov_temp);
for(InputIterator it = first;
it != beyond;
it++)
{
// Now for each triangle, construct the 2nd order moment about the origin.
// step 2: assemble the transformation matrix.
const Segment& t = *it;
// defined for convenience.
FT delta[4] = {t[0].x(), t[1].x(), t[0].y(), t[1].y()};
Matrix transformation = init_Matrix<K>(2,delta);
FT length = std::sqrt(t.squared_length());
CGAL_assertion(length!=0);
// step 3: Find the 2nd order moment for the triangle wrt to the origin by an affine transformation.
// step 3.1: Transform the standard 2nd order moment using the transformation matrix
transformation = length*transformation*covariance*LA::transpose(transformation);
//step 3.2: Translate the 2nd order moment to (x0,y0):Not Required Here
cov[0]+=transformation[0][0];
cov[1]+=transformation[0][1];
cov[2]+=transformation[1][1];
mass+=length;
}
//step 4: Translace the 2nd order miment calculated about the origin to the center of mass to get the covariance.
cov[0]+=mass*(-1.0*c.x()*c.x());
cov[1]+=mass*(-1.0*c.x()*c.y());
cov[2]+=mass*(-1.0*c.y()*c.y());
std::cout<<cov[0]<<" "<<cov[1]<<" "<<cov[2]<<std::endl;
// solve for eigenvalues and eigenvectors.
// eigen values are sorted in descending order,
// eigen vectors are sorted in accordance.
std::pair<FT,FT> eigen_values;
std::pair<Vector,Vector> eigen_vectors;
// CGALi::eigen_symmetric_2<K>(final_cov, eigen_vectors, eigen_values);
FT eigen_vectors1[4];
FT eigen_values1[2];
eigen_symmetric<FT>(cov,2, eigen_vectors1, eigen_values1);
eigen_values = std::make_pair(eigen_values1[0],eigen_values1[1]);
eigen_vectors = std::make_pair(Vector(eigen_vectors1[0],eigen_vectors1[1]),Vector(eigen_vectors1[2],eigen_vectors1[3]));
// check unicity and build fitting line accordingly
if(eigen_values.first != eigen_values.second)
{
// regular case
line = Line(c, eigen_vectors.first);
return (FT)1.0 - eigen_values.second / eigen_values.first;
}
else
{
// isotropic case (infinite number of directions)
// by default: assemble a line that goes through
// the centroid and with a default horizontal vector.
line = Line(c, Vector(1.0, 0.0));
return (FT)0.0;
}
} // end linear_least_squares_fitting_2 for triangle set
} // end namespace CGALi
CGAL_END_NAMESPACE
#endif

View File

@ -0,0 +1,154 @@
// Copyright (c) 2005 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you may redistribute it under
// the terms of the Q Public License version 1.0.
// See the file LICENSE.QPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL: svn+ssh://gankit@scm.gforge.inria.fr/svn/cgal/trunk/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles.h $
// $Id: linear_least_squares_fitting_2.h 37882 2007-04-03 15:15:30Z spion $
//
// Author(s) : Pierre Alliez and Sylvain Pion and Ankit Gupta
#ifndef CGAL_LINEAR_LEAST_SQUARES_FITTING_TRIANGLES_H
#define CGAL_LINEAR_LEAST_SQUARES_FITTING_TRIANGLES_H
#include <CGAL/basic.h>
#include <CGAL/Object.h>
#include <CGAL/centroid.h>
#include <CGAL/eigen_2.h>
#include <CGAL/eigen.h>
#include <CGAL/Linear_algebraCd.h>
#include <iterator>
#include <vector>
#include <cmath>
CGAL_BEGIN_NAMESPACE
namespace CGALi {
template<typename K>
typename CGAL::Linear_algebraCd<typename K::FT>::Matrix
init_Matrix(int n, typename K::FT entries[]) {
CGAL_assertion(n>1);
typedef typename CGAL::Linear_algebraCd<typename K::FT>::Matrix Matrix;
Matrix ret(n);
for(int i=0;i<n;i++) {
for(int j=0;j<n;j++) {
ret[i][j]=entries[i*n+j];
}
}
return ret;
}
// Fits a line to a 2D triangle set.
// Returns a fitting quality (1 - lambda_min/lambda_max):
// 1 is best (zero variance orthogonally to the fitting line);
// 0 is worst (isotropic case, returns a line with horizontal
// direction by default)
template < typename InputIterator, typename K >
typename K::FT
linear_least_squares_fitting_2(InputIterator first,
InputIterator beyond,
typename K::Line_2& line, // best fit line
typename K::Point_2& c, // centroid
const K&, // kernel
const typename K::Triangle_2*) // used for indirection
{
// types
typedef typename K::FT FT;
typedef typename K::Line_2 Line;
typedef typename K::Point_2 Point;
typedef typename K::Vector_2 Vector;
typedef typename K::Triangle_2 Triangle;
typedef typename CGAL::Linear_algebraCd<FT> LA;
typedef typename LA::Matrix Matrix;
// precondition: at least one element in the container.
CGAL_precondition(first != beyond);
// compute centroid
c = centroid(first,beyond,K());
// assemble covariance matrix as a semi-definite matrix.
// Matrix numbering:
// 0
// 1 2
//Final combined covariance matrix for all triangles and their combined mass
FT mass=0.0,cov[3]={0.0,0.0,0.0};
// step 1: assemble the 2nd order moment about the origin.
FT cov_temp[4] = {1/12.0,1/24.0,1/24.0,1/12.0};
Matrix covariance = init_Matrix<K>(2,cov_temp);
for(InputIterator it = first;
it != beyond;
it++)
{
// Now for each triangle, construct the 2nd order moment about the origin.
// step 2: assemble the transformation matrix.
const Triangle& t = *it;
// defined for convenience.
FT x0 = t[0].x(), y0 = t[0].y();
FT delta[4] = {t[1].x() - x0,t[2].x() - x0,t[1].y() - y0,t[2].y() - y0};
Matrix transformation = init_Matrix<K>(2,delta);
FT detA = fabs(LA::determinant(transformation));
CGAL_assertion(detA!=0);
// step 3: Find the 2nd order moment for the triangle wrt to the origin by an affine transformation.
// step 3.1: Transform the standard 2nd order moment using the transformation matrix
transformation = detA*transformation*covariance*LA::transpose(transformation);
//step 3.2: Translate the 2nd order moment to (x0,y0).
FT xav0 = (delta[0]+delta[1])/3.0, yav0 = (delta[2]+delta[3])/3.0;
cov[0]+=transformation[0][0]+0.5*detA*(x0*xav0*2+x0*x0);
cov[1]+=transformation[0][1]+0.5*detA*(x0*yav0+xav0*y0+x0*y0);
cov[2]+=transformation[1][1]+0.5*detA*(y0*yav0*2+y0*y0);
mass+=0.5*detA;
}
//step 4: Translace the 2nd order miment calculated about the origin to the center of mass to get the covariance.
cov[0]+=mass*(-1.0*c.x()*c.x());
cov[1]+=mass*(-1.0*c.x()*c.y());
cov[2]+=mass*(-1.0*c.y()*c.y());
// solve for eigenvalues and eigenvectors.
// eigen values are sorted in descending order,
// eigen vectors are sorted in accordance.
std::pair<FT,FT> eigen_values;
std::pair<Vector,Vector> eigen_vectors;
// CGALi::eigen_symmetric_2<K>(final_cov, eigen_vectors, eigen_values);
FT eigen_vectors1[4];
FT eigen_values1[2];
eigen_symmetric<FT>(cov,2, eigen_vectors1, eigen_values1);
eigen_values = std::make_pair(eigen_values1[0],eigen_values1[1]);
eigen_vectors = std::make_pair(Vector(eigen_vectors1[0],eigen_vectors1[1]),Vector(eigen_vectors1[2],eigen_vectors1[3]));
// check unicity and build fitting line accordingly
if(eigen_values.first != eigen_values.second)
{
// regular case
line = Line(c, eigen_vectors.first);
return (FT)1.0 - eigen_values.second / eigen_values.first;
}
else
{
// isotropic case (infinite number of directions)
// by default: assemble a line that goes through
// the centroid and with a default horizontal vector.
line = Line(c, Vector(1.0, 0.0));
return (FT)0.0;
}
} // end linear_least_squares_fitting_2 for triangle set
} // end namespace CGALi
CGAL_END_NAMESPACE
#endif