diff --git a/.gitattributes b/.gitattributes
index ef025dcddd1..a922867e445 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -1896,15 +1896,30 @@ 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.png -text svneol=unset#image/png
 Principal_component_analysis/doc_tex/Principal_component_analysis/teaserLeastSquaresFitting.png -text
-Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_circles.cpp -text
-Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_rectangles.cpp -text
-Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_segments.cpp -text
-Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_triangles.cpp -text
+Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_ankit.cpp -text
+Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_circles_2.cpp -text
+Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_cuboids_3.cpp -text
+Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_points_2.cpp -text
+Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_points_3.cpp -text
+Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_rectangles_2.cpp -text
+Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_segments_2.cpp -text
+Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_segments_3.cpp -text
+Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_spheres_3.cpp -text
+Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_tetrahedrons_3.cpp -text
+Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_triangles_2.cpp -text
+Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_triangles_3.cpp -text
+Principal_component_analysis/include/CGAL/PCA_tags.h -text
 Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h -text
+Principal_component_analysis/include/CGAL/linear_least_squares_fitting_cuboids_3.h -text
 Principal_component_analysis/include/CGAL/linear_least_squares_fitting_points_2.h -text
+Principal_component_analysis/include/CGAL/linear_least_squares_fitting_points_3.h -text
 Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h -text
 Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h -text
+Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_3.h -text
+Principal_component_analysis/include/CGAL/linear_least_squares_fitting_spheres_3.h -text
+Principal_component_analysis/include/CGAL/linear_least_squares_fitting_tetrahedrons_3.h -text
 Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h -text
+Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h -text
 Principal_component_analysis/include/CGAL/util.h -text
 Principal_component_analysis/test/Principal_component_analysis/fitting_segments_2.cpp -text
 QP_solver/doc_tex/QP_solver/PkgDescription.tex -text
diff --git a/Principal_component_analysis/demo/Principal_component_analysis/windows/3d/MeshDoc.cpp b/Principal_component_analysis/demo/Principal_component_analysis/windows/3d/MeshDoc.cpp
index 87cbe0b56b5..e2e9e111f7a 100644
--- a/Principal_component_analysis/demo/Principal_component_analysis/windows/3d/MeshDoc.cpp
+++ b/Principal_component_analysis/demo/Principal_component_analysis/windows/3d/MeshDoc.cpp
@@ -85,41 +85,41 @@ BOOL CMeshDoc::OnOpenDocument(LPCTSTR lpszPathName)
 	// off extension
 	if(extension == ".off")
 	{
-		// read from stream
-		std::ifstream stream(lpszPathName);
-		if(!stream)
-		{
-			AfxMessageBox("Unable to open file");
-			return false;
-		}
-		stream >> m_mesh;
-
-    // add mesh points to point set
-    for(Mesh::Point_iterator it = m_mesh.points_begin();
-        it != m_mesh.points_end();
-        it++)
-      m_points.push_back(*it);
-
-    // add mesh triangles to triangle set
-    for(Mesh::Facet_iterator f = m_mesh.facets_begin();
-        f != m_mesh.facets_end();
-        f++)
-    {
-      const Point& a = f->halfedge()->vertex()->point();
-      const Point& b = f->halfedge()->next()->vertex()->point();
-      const Point& c = f->halfedge()->next()->next()->vertex()->point();
-      m_triangles.push_back(Triangle(a,b,c));
-    }
+	  // read from stream
+	  std::ifstream stream(lpszPathName);
+	  if(!stream)
+	    {
+	      AfxMessageBox("Unable to open file");
+	      return false;
+	    }
+	  stream >> m_mesh;
+	  
+	  // add mesh points to point set
+	  for(Mesh::Point_iterator it = m_mesh.points_begin();
+	      it != m_mesh.points_end();
+	      it++)
+	    m_points.push_back(*it);
+	  
+	  // add mesh triangles to triangle set
+	  for(Mesh::Facet_iterator f = m_mesh.facets_begin();
+	      f != m_mesh.facets_end();
+	      f++)
+	    {
+	      const Point& a = f->halfedge()->vertex()->point();
+	      const Point& b = f->halfedge()->next()->vertex()->point();
+	      const Point& c = f->halfedge()->next()->next()->vertex()->point();
+	      m_triangles.push_back(Triangle(a,b,c));
+	    }
 	}
 	else
-		{
-			AfxMessageBox("Unknown extension");
-			return false;
-		}
-
-  OnFitFitpointset();
-
+	  {
+	    AfxMessageBox("Unknown extension");
+	    return false;
+	  }
+	
+	OnFitFitpointset();
 
+	
 	return TRUE;
 }
 
diff --git a/Principal_component_analysis/examples/Principal_component_analysis/centroid.cpp b/Principal_component_analysis/examples/Principal_component_analysis/centroid.cpp
index f09d7b805cc..458f830d2de 100644
--- a/Principal_component_analysis/examples/Principal_component_analysis/centroid.cpp
+++ b/Principal_component_analysis/examples/Principal_component_analysis/centroid.cpp
@@ -19,7 +19,7 @@ int main()
   points_2.push_back(Point_2(2.0, 2.0));
   points_2.push_back(Point_2(3.0, 5.0));
 
-  Point_2 c2 = CGAL::centroid(points_2.begin(), points_2.end());
+  Point_2 c2 = CGAL::centroid(points_2.begin(), points_2.end(),CGAL::PCA_dimension_0_tag());
   std::cout << c2 << std::endl;
 
   // centroid of 3D points
@@ -28,7 +28,7 @@ int main()
   points_3.push_back(Point_3(2.0, 2.0, 1.2));
   points_3.push_back(Point_3(3.0, 5.0, 4.5));
 
-  Point_3 c3 = CGAL::centroid(points_3.begin(), points_3.end());
+  Point_3 c3 = CGAL::centroid(points_3.begin(), points_3.end(),CGAL::PCA_dimension_0_tag());
   std::cout << c3 << std::endl;
 
   return 0;
diff --git a/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_ankit.cpp b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_ankit.cpp
new file mode 100644
index 00000000000..5156f4f0c08
--- /dev/null
+++ b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_ankit.cpp
@@ -0,0 +1,27 @@
+// Example program for the linear_least_square_fitting function
+
+#include <CGAL/Cartesian.h>
+#include <CGAL/linear_least_squares_fitting_2.h>
+
+#include <list>
+typedef double               FT;
+typedef CGAL::Cartesian<FT>  K;
+typedef K::Line_2            Line_2;
+typedef K::Point_2           Point_2;
+
+int main(int argc, char** argv)
+{
+  std::list<Point_2> points;
+  points.push_back(Point_2(1.0,0.0));
+  points.push_back(Point_2(0.0,1.0));
+  points.push_back(Point_2(1.0,1.0));
+  points.push_back(Point_2(-1.0,1.0));
+  points.push_back(Point_2(-1.0,0.0));
+  //  points.push_back(Point_2(0.0,1.0));
+
+  Line_2 line;
+  linear_least_squares_fitting_2(points.begin(),points.end(),line,CGAL::PCA_dimension_0_tag());
+  
+  std::cout<<line<<std::endl;
+  return 0;
+}
diff --git a/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_circles.cpp b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_circles_2.cpp
similarity index 80%
rename from Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_circles.cpp
rename to Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_circles_2.cpp
index af610eec5cd..76566f49d5f 100644
--- a/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_circles.cpp
+++ b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_circles_2.cpp
@@ -15,10 +15,11 @@ int main()
 {
   std::list<Circle_2> circles;
   circles.push_back(Circle_2(Point_2(0.0,0.0),9));
-  circles.push_back(Circle_2(Point_2(1.0,1.0),25));
+  circles.push_back(Circle_2(Point_2(0.0,10.0),49));
+  circles.push_back(Circle_2(Point_2(10.0,0.0),49));
 
   Line_2 line;
-  FT i = linear_least_squares_fitting_2(circles.begin(),circles.end(),line);
+  FT i = linear_least_squares_fitting_2(circles.begin(),circles.end(),line,CGAL::PCA_dimension_1_tag());
 
   std::cout<<"accuracy: "<<i<<std::endl;
   std::cout<<line<<std::endl;
diff --git a/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_cuboids_3.cpp b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_cuboids_3.cpp
new file mode 100644
index 00000000000..9baa084e3cc
--- /dev/null
+++ b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_cuboids_3.cpp
@@ -0,0 +1,35 @@
+// Example program for the linear_least_square_fitting function
+
+#include <CGAL/Cartesian.h>
+#include <CGAL/linear_least_squares_fitting_3.h>
+
+#include <list>
+
+typedef double               FT;
+typedef CGAL::Cartesian<FT>  K;
+typedef K::Line_3            Line_3;
+typedef K::Plane_3          Plane_3;
+typedef K::Point_3           Point_3;
+typedef K::Iso_cuboid_3   Iso_cuboid_3;
+
+int main()
+{
+  std::list<Iso_cuboid_3> Iso_cuboids;
+  Iso_cuboids.push_back(Iso_cuboid_3(Point_3(0.0,0.0,0.0),Point_3(1.0,1.0,1.0)));
+  //  Iso_cuboids.push_back(Iso_cuboid_3(Point_3(4.0,0.0),Point_3(5.0,1.0)));
+
+  Line_3 line;
+  FT i = linear_least_squares_fitting_3(Iso_cuboids.begin(),Iso_cuboids.end(),line,CGAL::PCA_dimension_3_tag());
+
+  std::cout<<"Line's accuracy: "<<i<<std::endl;
+  std::cout<<"Line direction: "<<line.to_vector()<<std::endl;
+  std::cout<<"Line's point: "<<line.point(0)<<std::endl;
+
+  Plane_3 plane;
+  FT j = linear_least_squares_fitting_3(Iso_cuboids.begin(),Iso_cuboids.end(),plane,CGAL::PCA_dimension_3_tag());
+  
+  std::cout<<"Plane's accuracy: "<<j<<std::endl;
+  std::cout<<"Plane: "<<plane<<std::endl;
+  
+  return 0;
+}
diff --git a/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting.cpp b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_points_2.cpp
similarity index 94%
rename from Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting.cpp
rename to Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_points_2.cpp
index 61194c0193f..4ec63f605d9 100644
--- a/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting.cpp
+++ b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_points_2.cpp
@@ -18,7 +18,7 @@ int main()
   points.push_back(Point_2(3.0,0.0));
 
   Line_2 line;
-  linear_least_squares_fitting_2(points.begin(),points.end(),line);
+  linear_least_squares_fitting_2(points.begin(),points.end(),line,CGAL::PCA_dimension_0_tag());
 
   return 0;
 }
diff --git a/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_points_3.cpp b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_points_3.cpp
new file mode 100644
index 00000000000..5f4331b15b3
--- /dev/null
+++ b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_points_3.cpp
@@ -0,0 +1,24 @@
+// Example program for the linear_least_square_fitting function
+
+#include <CGAL/Cartesian.h>
+#include <CGAL/linear_least_squares_fitting_3.h>
+
+#include <list>
+
+typedef double               FT;
+typedef CGAL::Cartesian<FT>  K;
+typedef K::Line_3            Line_3;
+typedef K::Point_3           Point_3;
+
+int main()
+{
+  std::list<Point_3> points;
+  points.push_back(Point_3(1.0,0.0,0.0));
+  points.push_back(Point_3(2.0,0.0,0.0));
+  points.push_back(Point_3(3.0,0.0,0.0));
+
+  Line_3 line;
+  linear_least_squares_fitting_3(points.begin(),points.end(),line,CGAL::PCA_dimension_0_tag());
+
+  return 0;
+}
diff --git a/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_rectangles.cpp b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_rectangles_2.cpp
similarity index 90%
rename from Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_rectangles.cpp
rename to Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_rectangles_2.cpp
index 50a2e6e1ee7..6ac225ed143 100644
--- a/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_rectangles.cpp
+++ b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_rectangles_2.cpp
@@ -14,11 +14,11 @@ typedef K::Iso_rectangle_2   Iso_rectangle_2;
 int main()
 {
   std::list<Iso_rectangle_2> Iso_rectangles;
-  Iso_rectangles.push_back(Iso_rectangle_2(Point_2(0.0,0.0),Point_2(1.0,1.0)));
+  Iso_rectangles.push_back(Iso_rectangle_2(Point_2(0.0,0.0),Point_2(4.0,4.0)));
   Iso_rectangles.push_back(Iso_rectangle_2(Point_2(4.0,0.0),Point_2(5.0,1.0)));
 
   Line_2 line;
-  FT i = linear_least_squares_fitting_2(Iso_rectangles.begin(),Iso_rectangles.end(),line);
+  FT i = linear_least_squares_fitting_2(Iso_rectangles.begin(),Iso_rectangles.end(),line,CGAL::PCA_dimension_0_tag());
 
   std::cout<<"accuracy: "<<i<<std::endl;
   std::cout<<line<<std::endl;
diff --git a/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_segments.cpp b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_segments_2.cpp
similarity index 64%
rename from Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_segments.cpp
rename to Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_segments_2.cpp
index 18fcc7cb57c..538ddde326a 100644
--- a/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_segments.cpp
+++ b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_segments_2.cpp
@@ -14,13 +14,16 @@ typedef K::Segment_2         Segment_2;
 int main()
 {
   std::list<Segment_2> segments;
-  segments.push_back(Segment_2(Point_2(0.0,0.0),Point_2(0.0,1.0)));
-  segments.push_back(Segment_2(Point_2(0.0,2.0),Point_2(0.0,10.0)));
+  //  segments.push_back(Segment_2(Point_2(0.0,0.0),Point_2(0.0,1.0)));
+  segments.push_back(Segment_2(Point_2(0.0,1.0),Point_2(1.0,0.0)));
+  // segments.push_back(Segment_2(Point_2(0.0,0.0),Point_2(1.0,0.0)));
 
   Line_2 line;
-  FT i = linear_least_squares_fitting_2(segments.begin(),segments.end(),line);
+  Point_2 c;
+  FT i = linear_least_squares_fitting_2(segments.begin(),segments.end(),line,c,CGAL::PCA_dimension_1_tag());
 
   std::cout<<"accuracy: "<<i<<std::endl;
   std::cout<<line<<std::endl;
+  std::cout<<"centroid: "<<c<<std::endl;
   return 0;
 }
diff --git a/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_segments_3.cpp b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_segments_3.cpp
new file mode 100644
index 00000000000..5b61e253cb5
--- /dev/null
+++ b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_segments_3.cpp
@@ -0,0 +1,37 @@
+// Example program for the linear_least_square_fitting function
+
+#include <CGAL/Cartesian.h>
+#include <CGAL/linear_least_squares_fitting_3.h>
+
+#include <list>
+
+typedef double               FT;
+typedef CGAL::Cartesian<FT>  K;
+typedef K::Line_3            Line_3;
+typedef K::Plane_3           Plane_3;
+typedef K::Point_3           Point_3;
+typedef K::Segment_3         Segment_3;
+
+int main()
+{
+  std::list<Segment_3> segments;
+  segments.push_back(Segment_3(Point_3(1.0,1.0,1.0),Point_3(2.0,2.0,2.0)));
+  //  segments.push_back(Segment_3(Point_3(0.0,1.0),Point_3(1.0,0.0)));
+  //  segments.push_back(Segment_3(Point_3(0.0,0.0),Point_3(1.0,0.0)));
+
+  Line_3 line;
+  Point_3 c;
+  FT i = linear_least_squares_fitting_3(segments.begin(),segments.end(),line,c,CGAL::PCA_dimension_1_tag());
+
+  std::cout<<"accuracy: "<<i<<std::endl;
+  std::cout<<line<<std::endl;
+  std::cout<<"centroid: "<<c<<std::endl;
+
+  Plane_3 plane;
+  i = linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,CGAL::PCA_dimension_1_tag());
+
+  std::cout<<"accuracy: "<<i<<std::endl;
+  std::cout<<plane<<std::endl;
+  std::cout<<"centroid: "<<c<<std::endl;
+  return 0;
+}
diff --git a/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_spheres_3.cpp b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_spheres_3.cpp
new file mode 100644
index 00000000000..8b91a59f7bb
--- /dev/null
+++ b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_spheres_3.cpp
@@ -0,0 +1,33 @@
+// Example program for the linear_least_square_fitting function
+
+#include <CGAL/Cartesian.h>
+#include <CGAL/linear_least_squares_fitting_3.h>
+
+#include <list>
+
+typedef double               FT;
+typedef CGAL::Cartesian<FT>  K;
+typedef K::Line_3            Line_3;
+typedef K::Plane_3           Plane_3;
+typedef K::Point_3           Point_3;
+typedef K::Sphere_3          Sphere_3;
+
+int main()
+{
+  std::list<Sphere_3> spheres;
+  spheres.push_back(Sphere_3(Point_3(0.0,0.0,0.0),9));
+  spheres.push_back(Sphere_3(Point_3(0.0,10.0,0.0),25));
+
+  Line_3 line;
+  FT i = linear_least_squares_fitting_3(spheres.begin(),spheres.end(),line,CGAL::PCA_dimension_3_tag());
+
+  std::cout<<"Line's accuracy: "<<i<<std::endl;
+  std::cout<<"Line: "<<line<<std::endl;
+
+  Plane_3 plane;
+  FT j = linear_least_squares_fitting_3(spheres.begin(),spheres.end(),plane,CGAL::PCA_dimension_3_tag());
+
+  std::cout<<"Plane's accuracy: "<<j<<std::endl;
+  std::cout<<"Plane: "<<plane<<std::endl;
+  return 0;
+}
diff --git a/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_tetrahedrons_3.cpp b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_tetrahedrons_3.cpp
new file mode 100644
index 00000000000..50bee77160e
--- /dev/null
+++ b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_tetrahedrons_3.cpp
@@ -0,0 +1,37 @@
+// Example program for the linear_least_square_fitting function
+
+#include <CGAL/Cartesian.h>
+#include <CGAL/linear_least_squares_fitting_3.h>
+
+#include <list>
+
+typedef double               FT;
+typedef CGAL::Cartesian<FT>  K;
+typedef K::Line_3            Line_3;
+typedef K::Plane_3           Plane_3;
+typedef K::Point_3           Point_3;
+typedef K::Tetrahedron_3          Tetrahedron_3;
+
+int main()
+{
+  std::list<Tetrahedron_3> tetrahedrons;
+  tetrahedrons.push_back(Tetrahedron_3(Point_3(0.0,0.0,0.0),Point_3(1.0,0.0,0.0),Point_3(0.0,1.0,0.0),Point_3(0.0,0.0,1.0)));
+
+  Line_3 line;
+  FT i = linear_least_squares_fitting_3(tetrahedrons.begin(),tetrahedrons.end(),line,CGAL::PCA_dimension_3_tag());
+
+  std::cout<<"Line's accuracy: "<<i<<std::endl;
+  std::cout<<"Line direction: "<<line.to_vector()<<std::endl;
+  std::cout<<"Line's point: "<<line.point(0)<<std::endl;
+
+  if(line.has_on(Point_3(2/3.0,0,0)) && line.has_on(Point_3(0,2/3.0,0)))
+    std::cout<<"YES!!"<<std::endl;
+
+  Plane_3 plane;
+  FT j = linear_least_squares_fitting_3(tetrahedrons.begin(),tetrahedrons.end(),plane,CGAL::PCA_dimension_3_tag());
+  
+  std::cout<<"Plane's accuracy: "<<j<<std::endl;
+  std::cout<<"Plane: "<<plane<<std::endl;
+  
+  return 0;
+}
diff --git a/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_triangles.cpp b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_triangles_2.cpp
similarity index 66%
rename from Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_triangles.cpp
rename to Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_triangles_2.cpp
index 69e03b24119..15e6ee823d2 100644
--- a/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_triangles.cpp
+++ b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_triangles_2.cpp
@@ -14,13 +14,15 @@ typedef K::Triangle_2        Triangle_2;
 int main()
 {
   std::list<Triangle_2> triangles;
-  triangles.push_back(Triangle_2(Point_2(1.0,0.0),Point_2(0.0,1.0),Point_2(1.0,1.0)));
-  triangles.push_back(Triangle_2(Point_2(-1.0,0.0),Point_2(0.0,1.0),Point_2(-1.0,1.0)));
+  Point_2 c;
+  triangles.push_back(Triangle_2(Point_2(0.0,0.0),Point_2(0.0,1.0),Point_2(1.0,0.0)));
+  //  triangles.push_back(Triangle_2(Point_2(-1.0,0.0),Point_2(0.0,1.0),Point_2(-1.0,1.0)));
 
   Line_2 line;
-  FT i = linear_least_squares_fitting_2(triangles.begin(),triangles.end(),line);
+  FT i = linear_least_squares_fitting_2(triangles.begin(),triangles.end(),line,c,CGAL::PCA_dimension_0_tag());
 
   std::cout<<"accuracy: "<<i<<std::endl;
   std::cout<<line<<std::endl;
+  std::cout<<"centroid: "<<c<<std::endl;
   return 0;
 }
diff --git a/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_triangles_3.cpp b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_triangles_3.cpp
new file mode 100644
index 00000000000..8acbe5cfa5f
--- /dev/null
+++ b/Principal_component_analysis/examples/Principal_component_analysis/linear_least_squares_fitting_triangles_3.cpp
@@ -0,0 +1,37 @@
+// Example program for the linear_least_square_fitting function
+
+#include <CGAL/Cartesian.h>
+#include <CGAL/linear_least_squares_fitting_3.h>
+
+#include <list>
+
+typedef double               FT;
+typedef CGAL::Cartesian<FT>  K;
+typedef K::Line_3            Line_3;
+typedef K::Plane_3           Plane_3;
+typedef K::Point_3           Point_3;
+typedef K::Triangle_3          Triangle_3;
+
+int main()
+{
+  std::list<Triangle_3> triangles;
+  triangles.push_back(Triangle_3(Point_3(1.0,0.0,0.0),Point_3(0.0,1.0,0.0),Point_3(0.0,0.0,0.0)));
+
+  Line_3 line;
+  FT i = linear_least_squares_fitting_3(triangles.begin(),triangles.end(),line,CGAL::PCA_dimension_2_tag());
+
+  std::cout<<"Line's accuracy: "<<i<<std::endl;
+  std::cout<<"Line direction: "<<line.to_vector()<<std::endl;
+  std::cout<<"Line's point: "<<line.point(0)<<std::endl;
+
+  if(line.has_on(Point_3(2/3.0,0,0)) && line.has_on(Point_3(0,2/3.0,0)))
+    std::cout<<"YES!!"<<std::endl;
+
+  Plane_3 plane;
+  FT j = linear_least_squares_fitting_3(triangles.begin(),triangles.end(),plane,CGAL::PCA_dimension_2_tag());
+  
+  std::cout<<"Plane's accuracy: "<<j<<std::endl;
+  std::cout<<"Plane: "<<plane<<std::endl;
+  
+  return 0;
+}
diff --git a/Principal_component_analysis/include/CGAL/PCA_tags.h b/Principal_component_analysis/include/CGAL/PCA_tags.h
new file mode 100644
index 00000000000..d7ac35e7791
--- /dev/null
+++ b/Principal_component_analysis/include/CGAL/PCA_tags.h
@@ -0,0 +1,47 @@
+// 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/tags.h $
+// $Id: tags.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_TAGS_H
+#define CGAL_LINEAR_LEAST_SQUARES_FITTING_TAGS_H
+
+CGAL_BEGIN_NAMESPACE
+
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//::These Tags are meant to represent the type of
+//::geometry we are dealing with. The CGAL kernel
+//::provides only for circles, triangles, spheres,
+//: rectangles and so on. But these can be solid,
+//::hollow, just outlined etc. These tags resolve
+//::between these different types of the same geometry
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+// Fully filled geometry object, such as disks in 2d and balls in 3d
+struct PCA_dimension_3_tag {};
+
+// Empty geometry objects, such as triangles in 2d and hollow cubes with just the planar facets in 3d
+struct PCA_dimension_2_tag {};
+
+// Only relevant for 3d objects, outlines of 3d geometry objects formed by just lines joing appropriate vetrices, such as outline of a cuboid.
+struct PCA_dimension_1_tag {};
+
+// For the vertices of objects.
+struct PCA_dimension_0_tag {};
+
+CGAL_END_NAMESPACE
+
+#endif // CGAL_LINEAR_LEAST_SQUARES_FITTING_TAGS_H
diff --git a/Principal_component_analysis/include/CGAL/centroid.h b/Principal_component_analysis/include/CGAL/centroid.h
index e2586df42d4..f9780352818 100644
--- a/Principal_component_analysis/include/CGAL/centroid.h
+++ b/Principal_component_analysis/include/CGAL/centroid.h
@@ -25,7 +25,9 @@
 #include <iterator>
 #include <CGAL/Kernel_traits.h>
 #include <CGAL/Kernel/Dimension_utils.h>
+#include <CGAL/PCA_tags.h>
 
+#include <list>
 // Functions to compute the centroid of N points.
 // Works in 2D and 3D.
 
@@ -38,6 +40,8 @@ CGAL_BEGIN_NAMESPACE
 
 namespace CGALi {
 
+//:::::::::: 2D Objects :::::::::::::::::::
+
 // computes the centroid of a 2D point set
 // takes an iterator range over K::Point_2
 template < typename InputIterator, 
@@ -46,7 +50,8 @@ typename K::Point_2
 centroid(InputIterator begin, 
          InputIterator end, 
          const K&,
-         const typename K::Point_2*)
+         const typename K::Point_2*,
+	 const CGAL::PCA_dimension_0_tag& tag)
 {
   typedef typename K::Vector_2 Vector;
   typedef typename K::FT FT;
@@ -63,40 +68,45 @@ centroid(InputIterator begin,
   return ORIGIN + v / (FT)nb_pts;
 }// end centroid of a 2D point set
 
-// computes the centroid of a 3D point set
-// takes an iterator range over K::Point_3
-template < typename InputIterator, 
-           typename K >
-typename K::Point_3
-centroid(InputIterator begin, 
-         InputIterator end, 
-         const K& ,
-         const typename K::Point_3*)
-{
-  typedef typename K::Vector_3 Vector;
-  typedef typename K::FT FT;
-
-  CGAL_precondition(begin != end);
-
-  Vector v = NULL_VECTOR;
-  unsigned int nb_pts = 0;
-  while (begin != end) 
-  {
-    v = v + (*begin++ - ORIGIN);
-    nb_pts++;
-  }
-  return ORIGIN + v / (FT)nb_pts;
-}// end centroid of a 3D point set
-
 // computes the centroid of a 2D segment set
 // takes an iterator range over K::Segment_2
+
+// centriod for 2D segment set with 0D tag
+template < typename InputIterator, 
+           typename K >
+typename K::Point_2
+centroid(InputIterator begin, 
+         InputIterator end, 
+         const K& k,
+         const typename K::Segment_2*,
+	 const CGAL::PCA_dimension_0_tag& tag)
+{
+  typedef typename K::Point_2  Point;
+  typedef typename K::Segment_2 Segment;
+  
+  CGAL_precondition(begin != end);
+  
+  std::list<Point> points;  
+  for(InputIterator it = begin;
+      it != end;
+      it++)
+  {
+    const Segment& s = *it;
+    points.push_back(s[0]);
+    points.push_back(s[1]);
+  } 
+  return centroid(points.begin(),points.end(),k,(Point*)NULL,tag);
+}// end centriod for 2D segment set with 0D tag
+
+// centriod for 2D segment set with 1D tag
 template < typename InputIterator, 
            typename K >
 typename K::Point_2
 centroid(InputIterator begin, 
          InputIterator end, 
          const K& ,
-         const typename K::Segment_2*)
+         const typename K::Segment_2*,
+	 const CGAL::PCA_dimension_1_tag& tag)
 {
   typedef typename K::FT       FT;
   typedef typename K::Vector_2 Vector;
@@ -119,17 +129,78 @@ centroid(InputIterator begin,
   }
   CGAL_assertion(sum_lengths != 0.0);
   return ORIGIN + v / sum_lengths;
-} // end centroid of a 2D segment set
+} // end centroid of a 2D segment set with 1D tag
 
 // computes the centroid of a 2D triangle set
 // takes an iterator range over K::Triangle_2
+
+// centriod for 2D triangle set with 0D tag
+template < typename InputIterator, 
+           typename K >
+typename K::Point_2
+centroid(InputIterator begin, 
+         InputIterator end, 
+         const K& k,
+         const typename K::Triangle_2*,
+	 const CGAL::PCA_dimension_0_tag& tag)
+{
+  typedef typename K::Triangle_2 Triangle;
+  typedef typename K::Point_2 Point;
+
+  CGAL_precondition(begin != end);
+  
+  std::list<Point> points;
+  for(InputIterator it = begin;
+      it != end;
+      it++)
+  {
+    const Triangle& triangle = *it;
+    points.push_back(triangle[0]);
+    points.push_back(triangle[1]);
+    points.push_back(triangle[2]);
+  }
+  return centroid(points.begin(),points.end(),k,(Point*)NULL,tag);
+
+} // end centroid of a 2D triangle set with 0D tag
+
+// centriod for 2D triangle set with 1D tag
+template < typename InputIterator, 
+           typename K >
+typename K::Point_2
+centroid(InputIterator begin, 
+         InputIterator end, 
+         const K& k,
+         const typename K::Triangle_2*,
+	 const CGAL::PCA_dimension_1_tag& tag)
+{
+  typedef typename K::Triangle_2 Triangle;
+  typedef typename K::Segment_2 Segment;
+
+  CGAL_precondition(begin != end);
+  
+  std::list<Segment> segments;
+  for(InputIterator it = begin;
+      it != end;
+      it++)
+  {
+    const Triangle& triangle = *it;
+    segments.push_back(triangle[0],triangle[1]);
+    segments.push_back(triangle[1],triangle[2]);
+    segments.push_back(triangle[2],triangle[0]);
+  }
+  return centroid(segments.begin(),segments.end(),k,(Segment*)NULL,tag);
+
+} // end centroid of a 2D triangle set with 1D tag
+
+// centriod for 2D triangle set with 2D tag
 template < typename InputIterator, 
            typename K >
 typename K::Point_2
 centroid(InputIterator begin, 
          InputIterator end, 
          const K& ,
-         const typename K::Triangle_2*)
+         const typename K::Triangle_2*,
+	 const CGAL::PCA_dimension_2_tag& tag)
 {
   typedef typename K::FT       FT;
   typedef typename K::Vector_2 Vector;
@@ -152,17 +223,53 @@ centroid(InputIterator begin,
   }
   CGAL_assertion(sum_areas != 0.0);
   return ORIGIN + v / sum_areas;
-} // end centroid of a 2D triangle set
+} // end centroid of a 2D triangle set with 2D tag
 
 // computes the centroid of a 2D circle set
 // takes an iterator range over K::Circle_2
+
+// centriod for 2D circle set with 1D tag
 template < typename InputIterator, 
            typename K >
 typename K::Point_2
 centroid(InputIterator begin, 
          InputIterator end, 
          const K& ,
-         const typename K::Circle_2*)
+         const typename K::Circle_2*,
+	 const CGAL::PCA_dimension_1_tag& tag)
+{
+  typedef typename K::FT       FT;
+  typedef typename K::Vector_2 Vector;
+  typedef typename K::Point_2  Point;
+  typedef typename K::Circle_2 Circle;
+
+  CGAL_precondition(begin != end);
+
+  Vector v = NULL_VECTOR;
+  FT sum_lengths = 0;
+  for(InputIterator it = begin;
+      it != end;
+      it++)
+  {
+    const Circle& s = *it;
+    FT radius = std::sqrt(s.squared_radius());
+    Point c = s.center();
+    v = v + radius * (c - ORIGIN);
+    sum_lengths += radius;
+  }
+  CGAL_assertion(sum_lengths != 0.0);
+  return ORIGIN + v / sum_lengths;
+} // end centroid of a 2D circle set with 1D tag
+
+// centriod for 2D circle set with 2D tag
+template < typename InputIterator, 
+           typename K >
+typename K::Point_2
+centroid(InputIterator begin, 
+         InputIterator end, 
+         const K& ,
+         const typename K::Circle_2*,
+	 const CGAL::PCA_dimension_2_tag& tag)
 {
   typedef typename K::FT       FT;
   typedef typename K::Vector_2 Vector;
@@ -179,24 +286,86 @@ centroid(InputIterator begin,
   {
     const Circle& s = *it;
     FT sq_radius = s.squared_radius();
-    //    Point c = K().construct_midpoint_2_object()(s[0],s[1]);
     Point c = s.center();
     v = v + sq_radius * (c - ORIGIN);
     sum_areas += sq_radius;
   }
   CGAL_assertion(sum_areas != 0.0);
   return ORIGIN + v / sum_areas;
-} // end centroid of a 2D circle set
+} // end centroid of a 2D circle set with 2D tag
 
 // computes the centroid of a 2D rectangle set
 // takes an iterator range over K::Iso_Rectangle_2
+
+// centriod for 2D rectangle set with 0D tag
+template < typename InputIterator, 
+           typename K >
+typename K::Point_2
+centroid(InputIterator begin, 
+         InputIterator end, 
+         const K& k,
+         const typename K::Iso_rectangle_2*,
+	 const CGAL::PCA_dimension_0_tag& tag)
+{
+  typedef typename K::Iso_rectangle_2 Iso_rectangle;
+  typedef typename K::Point_2 Point;
+
+  CGAL_precondition(begin != end);
+  
+  std::list<Point> points;
+  for(InputIterator it = begin;
+      it != end;
+      it++)
+  {
+    const Iso_rectangle& r = *it;
+    points.push_back(r[0]);
+    points.push_back(r[1]);
+    points.push_back(r[2]);
+    points.push_back(r[3]);
+  }
+  return centroid(points.begin(),points.end(),k,(Point*)NULL,tag);
+
+} // end centroid of a 2D rectangle set with 0D tag
+
+// centriod for 2D rectangle set with 1D tag
+template < typename InputIterator, 
+           typename K >
+typename K::Point_2
+centroid(InputIterator begin, 
+         InputIterator end, 
+         const K& k,
+         const typename K::Iso_rectangle_2*,
+	 const CGAL::PCA_dimension_1_tag& tag)
+{
+  typedef typename K::Iso_rectangle_2 Iso_rectangle;
+  typedef typename K::Segment_2 Segment;
+
+  CGAL_precondition(begin != end);
+  
+  std::list<Segment> segments;
+  for(InputIterator it = begin;
+      it != end;
+      it++)
+  {
+    const Iso_rectangle& r = *it;
+    segments.push_back(r[0],r[1]);
+    segments.push_back(r[1],r[2]);
+    segments.push_back(r[2],r[3]);
+    segments.push_back(r[3],r[0]);
+  }
+  return centroid(segments.begin(),segments.end(),k,(Segment*)NULL,tag);
+
+} // end centroid of a 2D rectangle set with 1D tag
+
+// centriod for 2D rectangle set with 2D tag
 template < typename InputIterator, 
            typename K >
 typename K::Point_2
 centroid(InputIterator begin, 
          InputIterator end, 
          const K& ,
-         const typename K::Iso_rectangle_2*)
+         const typename K::Iso_rectangle_2*,
+	 const CGAL::PCA_dimension_2_tag& tag)
 {
   typedef typename K::FT       FT;
   typedef typename K::Vector_2 Vector;
@@ -219,17 +388,141 @@ centroid(InputIterator begin,
   }
   CGAL_assertion(sum_areas != 0.0);
   return ORIGIN + v / sum_areas;
-} // end centroid of a 2D rectangle set
+} // end centroid of a 2D rectangle set with 2D tag
 
-// computes the centroid of a 3D triangle set
-// takes an iterator range over K::Triangle_3
+//::::::::::: 3D Objects ::::::::::::::::::::::::::
+
+// computes the centroid of a 3D point set
+// takes an iterator range over K::Point_3
+
+// centriod for 3D point set with 0D tag
+template < typename InputIterator, 
+           typename K >
+typename K::Point_3
+centroid(InputIterator begin, 
+         InputIterator end, 
+         const K&,
+         const typename K::Point_3*,
+	 const CGAL::PCA_dimension_0_tag& tag)
+{
+  typedef typename K::Vector_3 Vector;
+  typedef typename K::FT FT;
+
+  CGAL_precondition(begin != end);
+
+  Vector v = NULL_VECTOR;
+  unsigned int nb_pts = 0;
+  while (begin != end) 
+  {
+    v = v + (*begin++ - ORIGIN);
+    nb_pts++;
+  }
+  return ORIGIN + v / (FT)nb_pts;
+}// end centroid of a 3D point set with 0D tag
+
+// centriod for 3D segment set with 1D tag
 template < typename InputIterator, 
            typename K >
 typename K::Point_3
 centroid(InputIterator begin, 
          InputIterator end, 
          const K& ,
-         const typename K::Triangle_3*)
+         const typename K::Segment_3*,
+	 const CGAL::PCA_dimension_1_tag& tag)
+{
+  typedef typename K::FT       FT;
+  typedef typename K::Vector_3 Vector;
+  typedef typename K::Point_3  Point;
+  typedef typename K::Segment_3 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_midpoint_3_object()(s[0],s[1]);
+    v = v + length * (c - ORIGIN);
+    sum_lengths += length;
+  }
+  CGAL_assertion(sum_lengths != 0.0);
+  return ORIGIN + v / sum_lengths;
+} // end centroid of a 3D segment set with 1D tag
+
+// computes the centroid of a 3D triangle set
+// takes an iterator range over K::Triangle_3
+
+// centriod for 3D triangle set with 0D tag
+template < typename InputIterator, 
+           typename K >
+typename K::Point_3
+centroid(InputIterator begin, 
+         InputIterator end, 
+         const K& k,
+         const typename K::Triangle_3*,
+	 const CGAL::PCA_dimension_0_tag& tag)
+{
+  typedef typename K::Triangle_3 Triangle;
+  typedef typename K::Point_3 Point;
+
+  CGAL_precondition(begin != end);
+  
+  std::list<Point> points;
+  for(InputIterator it = begin;
+      it != end;
+      it++)
+  {
+    const Triangle& triangle = *it;
+    points.push_back(triangle[0]);
+    points.push_back(triangle[1]);
+    points.push_back(triangle[2]);
+  }
+  return centroid(points.begin(),points.end(),k,(Point*)NULL,tag);
+
+} // end centroid of a 3D triangle set with 0D tag
+
+// centriod for 3D triangle set with 1D tag
+template < typename InputIterator, 
+           typename K >
+typename K::Point_3
+centroid(InputIterator begin, 
+         InputIterator end, 
+         const K& k,
+         const typename K::Triangle_3*,
+	 const CGAL::PCA_dimension_1_tag& tag)
+{
+  typedef typename K::Triangle_3 Triangle;
+  typedef typename K::Segment_3 Segment;
+
+  CGAL_precondition(begin != end);
+  
+  std::list<Segment> segments;
+  for(InputIterator it = begin;
+      it != end;
+      it++)
+  {
+    const Triangle& triangle = *it;
+    segments.push_back(triangle[0],triangle[1]);
+    segments.push_back(triangle[1],triangle[2]);
+    segments.push_back(triangle[2],triangle[0]);
+  }
+  return centroid(segments.begin(),segments.end(),k,(Segment*)NULL,tag);
+
+} // end centroid of a 3D triangle set with 1D tag
+
+// centriod for 3D triangle set with 2D tag
+template < typename InputIterator, 
+           typename K >
+typename K::Point_3
+centroid(InputIterator begin, 
+         InputIterator end, 
+         const K& ,
+         const typename K::Triangle_3*,
+	 const CGAL::PCA_dimension_2_tag& tag)
 {
   typedef typename K::FT       FT;
   typedef typename K::Vector_3 Vector;
@@ -252,39 +545,284 @@ centroid(InputIterator begin,
   }
   CGAL_assertion(sum_areas != 0.0);
   return ORIGIN + v / sum_areas;
-} // end centroid of a 3D triangle set
+} // end centroid of a 3D triangle set with 2D tag
+
+// computes the centroid of a 3D sphere set
+// takes an iterator range over K::Sphere_3
+
+// centriod for 3D sphere set with 2D tag
+template < typename InputIterator, 
+           typename K >
+typename K::Point_3
+centroid(InputIterator begin, 
+         InputIterator end, 
+         const K& ,
+         const typename K::Sphere_3*,
+	 const CGAL::PCA_dimension_2_tag& tag)
+{
+  typedef typename K::FT       FT;
+  typedef typename K::Vector_3 Vector;
+  typedef typename K::Point_3  Point;
+  typedef typename K::Sphere_3 Sphere;
+
+  CGAL_precondition(begin != end);
+
+  Vector v = NULL_VECTOR;
+  FT sum_areas = 0;
+  for(InputIterator it = begin;
+      it != end;
+      it++)
+  {
+    const Sphere& sphere = *it;
+    FT unsigned_area = sphere.squared_radius();
+    Point c = sphere.center();
+    v = v + unsigned_area * (c - ORIGIN);
+    sum_areas += unsigned_area;
+  }
+  CGAL_assertion(sum_areas != 0.0);
+  return ORIGIN + v / sum_areas;
+} // end centroid of a 3D sphere set with 2D tag
+
+// centriod for 3D sphere set with 3D tag
+template < typename InputIterator, 
+           typename K >
+typename K::Point_3
+centroid(InputIterator begin, 
+         InputIterator end, 
+         const K& ,
+         const typename K::Sphere_3*,
+	 const CGAL::PCA_dimension_3_tag& tag)
+{
+  typedef typename K::FT       FT;
+  typedef typename K::Vector_3 Vector;
+  typedef typename K::Point_3  Point;
+  typedef typename K::Sphere_3 Sphere;
+
+  CGAL_precondition(begin != end);
+
+  Vector v = NULL_VECTOR;
+  FT sum_volumes = 0;
+  for(InputIterator it = begin;
+      it != end;
+      it++)
+  {
+    const Sphere& sphere = *it;
+    FT unsigned_volume = sphere.squared_radius() * std::sqrt(sphere.squared_radius());
+    Point c = sphere.center();
+    v = v + unsigned_volume * (c - ORIGIN);
+    sum_volumes += unsigned_volume;
+  }
+  CGAL_assertion(sum_volumes != 0.0);
+  return ORIGIN + v / sum_volumes;
+} // end centroid of a 3D sphere set with 3 tag
+
+// computes the centroid of a 3D cuboid set
+// takes an iterator range over K::Iso_cuboid_3
+
+// centriod for 3D cuboid set with 0D tag
+template < typename InputIterator, 
+           typename K >
+typename K::Point_3
+centroid(InputIterator begin, 
+         InputIterator end, 
+          const K& k,
+         const typename K::Iso_cuboid_3*,
+	 const CGAL::PCA_dimension_0_tag& tag)
+{
+  typedef typename K::Iso_cuboid_3 Iso_cuboid;
+  typedef typename K::Point_3 Point;
+
+  CGAL_precondition(begin != end);
+  
+  std::list<Point> points;
+  for(InputIterator it = begin;
+      it != end;
+      it++)
+  {
+    const Iso_cuboid& cuboid = *it;
+    points.push_back(cuboid[0]);
+    points.push_back(cuboid[1]);
+    points.push_back(cuboid[2]);
+    points.push_back(cuboid[3]);
+    points.push_back(cuboid[4]);
+    points.push_back(cuboid[5]);
+    points.push_back(cuboid[6]);
+    points.push_back(cuboid[7]);
+  }
+  return centroid(points.begin(),points.end(),k,(Point*)NULL,tag);
+
+} // end centroid of a 3D cuboid set with 0D tag
+
+// centriod for 3D cuboid set with 1D tag
+template < typename InputIterator, 
+           typename K >
+typename K::Point_3
+centroid(InputIterator begin, 
+         InputIterator end, 
+          const K& k,
+         const typename K::Iso_cuboid_3*,
+	 const CGAL::PCA_dimension_1_tag& tag)
+{
+  typedef typename K::Iso_cuboid_3 Iso_cuboid;
+  typedef typename K::Segment_3 Segment;
+
+  CGAL_precondition(begin != end);
+  
+  std::list<Segment> segments;
+  for(InputIterator it = begin;
+      it != end;
+      it++)
+  {
+    const Iso_cuboid& cuboid = *it;
+    segments.push_back(cuboid[0],cuboid[1]);
+    segments.push_back(cuboid[1],cuboid[2]);
+    segments.push_back(cuboid[2],cuboid[3]);
+    segments.push_back(cuboid[3],cuboid[0]);
+    segments.push_back(cuboid[0],cuboid[5]);
+    segments.push_back(cuboid[5],cuboid[4]);
+    segments.push_back(cuboid[4],cuboid[3]);
+    segments.push_back(cuboid[1],cuboid[6]);
+    segments.push_back(cuboid[6],cuboid[7]);
+    segments.push_back(cuboid[7],cuboid[2]);
+    segments.push_back(cuboid[4],cuboid[7]);
+    segments.push_back(cuboid[5],cuboid[6]);
+  }
+  return centroid(segments.begin(),segments.end(),k,(Segment*)NULL,tag);
+
+} // end centroid of a 3D cuboid set with 1D tag
+
+// centriod for 3D cuboid set with 2D tag
+template < typename InputIterator, 
+           typename K >
+typename K::Point_3
+centroid(InputIterator begin, 
+         InputIterator end, 
+          const K& ,
+         const typename K::Iso_cuboid_3*,
+	 const CGAL::PCA_dimension_2_tag& tag)
+{
+  typedef typename K::FT       FT;
+  typedef typename K::Vector_3 Vector;
+  typedef typename K::Point_3  Point;
+  typedef typename K::Iso_cuboid_3 Iso_cuboid;
+
+  CGAL_precondition(begin != end);
+
+  Vector v = NULL_VECTOR;
+  FT sum_areas = 0;
+  for(InputIterator it = begin;
+      it != end;
+      it++)
+  {
+    const Iso_cuboid& cuboid = *it;
+    FT unsigned_area = 2 * ((cuboid.xmax()-cuboid.xmin())*(cuboid.ymax()-cuboid.ymin()) + (cuboid.xmax()-cuboid.xmin())*(cuboid.zmax()-cuboid.zmin()) + (cuboid.ymax()-cuboid.ymin())*(cuboid.zmax()-cuboid.zmin()));
+    Point c = K().construct_centroid_3_object()(cuboid[0],cuboid[1],cuboid[3],cuboid[5]);
+    v = v + unsigned_area * (c - ORIGIN);
+    sum_areas += unsigned_area;
+  }
+  CGAL_assertion(sum_areas != 0.0);
+  return ORIGIN + v / sum_areas;
+} // end centroid of a 3D cuboid set with 2D tag
+
+// centriod for 3D cuboid set with 3D tag
+template < typename InputIterator, 
+           typename K >
+typename K::Point_3
+centroid(InputIterator begin, 
+         InputIterator end, 
+          const K& ,
+         const typename K::Iso_cuboid_3*,
+	 const CGAL::PCA_dimension_3_tag& tag)
+{
+  typedef typename K::FT       FT;
+  typedef typename K::Vector_3 Vector;
+  typedef typename K::Point_3  Point;
+  typedef typename K::Iso_cuboid_3 Iso_cuboid;
+
+  CGAL_precondition(begin != end);
+
+  Vector v = NULL_VECTOR;
+  FT sum_volumes = 0;
+  for(InputIterator it = begin;
+      it != end;
+      it++)
+  {
+    const Iso_cuboid& cuboid = *it;
+    FT unsigned_volume = cuboid.volume();
+    Point c = K().construct_centroid_3_object()(cuboid[0],cuboid[1],cuboid[3],cuboid[5]);
+    v = v + unsigned_volume * (c - ORIGIN);
+    sum_volumes += unsigned_volume;
+  }
+  CGAL_assertion(sum_volumes != 0.0);
+  return ORIGIN + v / sum_volumes;
+} // end centroid of a 3D cuboid set with 3D tag
+
+// centriod for 3D Tetrahedron set with 3D tag
+template < typename InputIterator, 
+           typename K >
+typename K::Point_3
+centroid(InputIterator begin, 
+         InputIterator end, 
+          const K& ,
+         const typename K::Tetrahedron_3*,
+	 const CGAL::PCA_dimension_3_tag& tag)
+{
+  typedef typename K::FT       FT;
+  typedef typename K::Vector_3 Vector;
+  typedef typename K::Point_3  Point;
+  typedef typename K::Tetrahedron_3 Tetrahedron;
+
+  CGAL_precondition(begin != end);
+
+  Vector v = NULL_VECTOR;
+  FT sum_volumes = 0;
+  for(InputIterator it = begin;
+      it != end;
+      it++)
+  {
+    const Tetrahedron& Tetrahedron = *it;
+    FT unsigned_volume = Tetrahedron.volume();
+    Point c = K().construct_centroid_3_object()(Tetrahedron[0],Tetrahedron[1],Tetrahedron[2],Tetrahedron[3]);
+    v = v + unsigned_volume * (c - ORIGIN);
+    sum_volumes += unsigned_volume;
+  }
+  CGAL_assertion(sum_volumes != 0.0);
+  return ORIGIN + v / sum_volumes;
+} // end centroid of a 3D Tetrahedron set with 3D tag
 
 } // namespace CGALi
 
 // computes the centroid of a set of kernel objects
 // takes an iterator range over kernel objects
 template < typename InputIterator, 
-           typename K >
+           typename K, 
+	   typename tag  >
 inline
 typename Point<Dimension<typename std::iterator_traits<InputIterator>::value_type, K>::value,
                K
               >::type
 centroid(InputIterator begin,
          InputIterator end, 
-         const K& k)
+         const K& k,
+	 const tag& t)
 {
   typedef typename std::iterator_traits<InputIterator>::value_type Value_type;
-  return CGALi::centroid(begin, end, k,(Value_type*) NULL);
+  return CGALi::centroid(begin, end, k,(Value_type*) NULL, t);
 }
 
 // this one takes an iterator range over kernel objects
 // and uses Kernel_traits<> to find out its kernel.
-template < typename InputIterator >
+template < typename InputIterator, typename tag >
 inline
 typename Point<Dimension<
                typename std::iterator_traits<InputIterator>::value_type,
                typename Kernel_traits<typename std::iterator_traits<InputIterator>::value_type>::Kernel >::value,
                typename Kernel_traits<typename std::iterator_traits<InputIterator>::value_type>::Kernel >::type
-centroid(InputIterator begin, InputIterator end)
+centroid(InputIterator begin, InputIterator end, const tag& t)
 {
   typedef typename std::iterator_traits<InputIterator>::value_type  Point;
   typedef typename Kernel_traits<Point>::Kernel                     K;
-  return CGAL::centroid(begin, end, K());
+  return CGAL::centroid(begin, end, K(), t);
 }
 
 CGAL_END_NAMESPACE
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_2.h
index a788ab07783..54def822d72 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_2.h
@@ -21,19 +21,21 @@
 
 #include <CGAL/basic.h>
 #include <CGAL/Object.h>
-#include <CGAL/centroid.h>
+#include <CGAL/Algebraic_structure_traits.h>
+#include <CGAL/IO/io.h>
 #include <CGAL/linear_least_squares_fitting_points_2.h>
 #include <CGAL/linear_least_squares_fitting_segments_2.h>
 #include <CGAL/linear_least_squares_fitting_triangles_2.h>
 #include <CGAL/linear_least_squares_fitting_circles_2.h>
 #include <CGAL/linear_least_squares_fitting_rectangles_2.h>
+#include <CGAL/PCA_tags.h>
 
 #include <iterator>
 
 CGAL_BEGIN_NAMESPACE
 
 template < typename InputIterator, 
-           typename K >
+           typename K , typename tag>
 inline
 typename K::FT
 linear_least_squares_fitting_2(InputIterator first,
@@ -41,59 +43,63 @@ linear_least_squares_fitting_2(InputIterator first,
                                typename K::Line_2& line,
                                typename K::Point_2& centroid,
                                const K& k,
-			       const bool non_standard_geometry = false)
+			       const tag& t)
 {
   typedef typename std::iterator_traits<InputIterator>::value_type Value_type;
+  //  BOOST_STATIC_ASSERT((boost::is_same<typename CGAL::Algebraic_structure_traits<Value_type>::Algebraic_category,CGAL::Field_with_sqrt_tag>::value));
   return CGALi::linear_least_squares_fitting_2(first, beyond, line,
-                                               centroid, k, (Value_type*) NULL, non_standard_geometry);
+                                               centroid, k, (Value_type*) NULL, t);
 }
 
 template < typename InputIterator, 
-           typename K >
+           typename K, typename tag >
 inline
 typename K::FT
 linear_least_squares_fitting_2(InputIterator first,
                                InputIterator beyond, 
                                typename K::Line_2& line,
                                const K& k,
-			       bool non_standard_geometry = false)
+			       const tag& t)
 {
   typedef typename std::iterator_traits<InputIterator>::value_type Value_type;
+  //  BOOST_STATIC_ASSERT((boost::is_same<typename CGAL::Algebraic_structure_traits<Value_type>::Algebraic_category,CGAL::Field_with_sqrt_tag>::value));
   typename K::Point_2 centroid;
   return CGALi::linear_least_squares_fitting_2(first, beyond, line,
-                                               centroid, k,(Value_type*) NULL, non_standard_geometry);
+                                               centroid, k,(Value_type*) NULL, t);
 }
 
 // deduces the kernel from the points in container.
 template < typename InputIterator, 
            typename Line,
-           typename Point>
+           typename Point, typename tag>
 inline
 typename Kernel_traits<Line>::Kernel::FT
 linear_least_squares_fitting_2(InputIterator first,
                                InputIterator beyond, 
                                Line& line,
                                Point& centroid,
-			       bool non_standard_geometry = false)
+			       const tag& t)
 {
   typedef typename std::iterator_traits<InputIterator>::value_type Value_type;
+  //  BOOST_STATIC_ASSERT((boost::is_same<typename CGAL::Algebraic_structure_traits<Value_type>::Algebraic_category,CGAL::Field_with_sqrt_tag>::value));
   typedef typename Kernel_traits<Value_type>::Kernel K;
-  return CGAL::linear_least_squares_fitting_2(first,beyond,line,centroid,K(), non_standard_geometry);
+  return CGAL::linear_least_squares_fitting_2(first,beyond,line,centroid,K(), t);
 }
 
 // does not return the centroid and deduces the kernel as well.
 template < typename InputIterator, 
-           typename Line >
+           typename Line, typename tag >
 inline
 typename Kernel_traits<Line>::Kernel::FT
 linear_least_squares_fitting_2(InputIterator first,
                                InputIterator beyond, 
                                Line& line,
-			       bool non_standard_geometry = false)
+			       const tag& t)
 {
   typedef typename std::iterator_traits<InputIterator>::value_type Value_type;
   typedef typename Kernel_traits<Value_type>::Kernel K;
-  return CGAL::linear_least_squares_fitting_2(first,beyond,line,K(), non_standard_geometry);
+  //  BOOST_STATIC_ASSERT((boost::is_same<typename CGAL::Algebraic_structure_traits<Value_type>::Algebraic_category,CGAL::Field_with_sqrt_tag>::value));
+  return CGAL::linear_least_squares_fitting_2(first,beyond,line,K(), t);
 }
 
 CGAL_END_NAMESPACE
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_3.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_3.h
index 6dc7b93e062..1ce6e2a8de5 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_3.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_3.h
@@ -14,15 +14,22 @@
 // $URL$
 // $Id$
 //
-// Author(s)     : Pierre Alliez and Sylvain Pion
+// Author(s)     : Pierre Alliez and Sylvain Pion and Ankit Gupta
 
 #ifndef CGAL_LINEAR_LEAST_SQUARES_FITTING_3_H
 #define CGAL_LINEAR_LEAST_SQUARES_FITTING_3_H
 
 #include <CGAL/basic.h>
 #include <CGAL/Object.h>
-#include <CGAL/centroid.h>
-#include <CGAL/eigen.h>
+#include <CGAL/Algebraic_structure_traits.h>
+#include <CGAL/IO/io.h>
+#include <CGAL/linear_least_squares_fitting_points_3.h>
+#include <CGAL/linear_least_squares_fitting_segments_3.h>
+#include <CGAL/linear_least_squares_fitting_triangles_3.h>
+#include <CGAL/linear_least_squares_fitting_spheres_3.h>
+#include <CGAL/linear_least_squares_fitting_cuboids_3.h>
+#include <CGAL/linear_least_squares_fitting_tetrahedrons_3.h>
+#include <CGAL/PCA_tags.h>
 
 #include <iterator>
 #include <list>
@@ -30,413 +37,106 @@
 
 CGAL_BEGIN_NAMESPACE
 
-namespace CGALi {
-
-
-// assemble covariance matrix from a point set 
-template < typename InputIterator,
-           typename K >
-void
-assemble_covariance_matrix_3(InputIterator first,
-                             InputIterator beyond, 
-                             typename K::FT covariance[6], // covariance matrix
-                             const typename K::Point_3& c, // centroid
-                             const K& ,                    // kernel
-                             const typename K::Point_3*)   // used for indirection
-{
-  typedef typename K::FT       FT;
-  typedef typename K::Point_3  Point;
-  typedef typename K::Vector_3 Vector;
-
-  // Matrix numbering:
-  // 0          
-  // 1 2
-  // 3 4 5          
-  covariance[0] = covariance[1] = covariance[2] = 
-  covariance[3] = covariance[4] = covariance[5] = (FT)0.0;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Point& p = *it;
-    Vector d = p - c;
-    covariance[0] += d.x() * d.x();
-    covariance[1] += d.x() * d.y();
-    covariance[2] += d.y() * d.y();
-    covariance[3] += d.x() * d.z();
-    covariance[4] += d.y() * d.z();
-    covariance[5] += d.z() * d.z();
-  }
-}
-
-// assemble covariance matrix from a triangle set 
-template < typename InputIterator,
-           typename K >
-void
-assemble_covariance_matrix_3(InputIterator first,
-                             InputIterator beyond, 
-                             typename K::FT covariance[6], // covariance matrix
-                             const typename K::Point_3& c, // centroid
-                             const K& ,                    // kernel
-                             const typename K::Triangle_3*)// used for indirection
-{
-  typedef typename K::FT       FT;
-  typedef typename K::Point_3  Point;
-  typedef typename K::Vector_3 Vector;
-  typedef typename K::Triangle_3 Triangle;
-
-  // Matrix numbering:
-  // 0          
-  // 1 2
-  // 3 4 5          
-  covariance[0] = covariance[1] = covariance[2] = 
-  covariance[3] = covariance[4] = covariance[5] = (FT)0.0;
-  FT sum_areas = 0.0;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Triangle& triangle = *it;
-    FT area = std::sqrt(triangle.squared_area());
-    Point c_t = K().construct_centroid_3_object()(triangle[0],triangle[1],triangle[2]);
-    sum_areas += area;
-
-    // e1 = ab, e2 = ac
-    Vector e1 = Vector(triangle[0],triangle[1]);
-    Vector e2 = Vector(triangle[0],triangle[2]);
-
-    FT c1 = (FT)(2.0 * area * 10.0 / 72.0);
-    FT c2 = (FT)(2.0 * area *  7.0 / 72.0);
-        
-    covariance[0] += c1*(e1[0]*e1[0] + e2[0]*e2[0]) + (FT)2.0*c2*e1[0]*e2[0];
-    covariance[1] += c1*(e1[1]*e1[0] + e2[1]*e2[0]) + c2*(e1[1]*e2[0] + e1[0]*e2[1]);
-    covariance[2] += c1*(e1[1]*e1[1] + e2[1]*e2[1]) + (FT)2.0*c2*e1[1]*e2[1];
-    covariance[3] += c1*(e1[2]*e1[0] + e2[2]*e2[0]) + c2*(e1[2]*e2[0] + e1[0]*e2[2]);
-    covariance[4] += c1*(e1[2]*e1[1] + e2[2]*e2[1]) + c2*(e1[2]*e2[1] + e1[1]*e2[2]);
-    covariance[5] += c1*(e1[2]*e1[2] + e2[2]*e2[2]) + (FT)2.0*c2*e1[2]*e2[2];
-    
-    // add area(t) c(t) * transpose(c(t))
-    covariance[0] += area * c_t.x() * c_t.x();
-    covariance[1] += area * c_t.y() * c_t.x();
-    covariance[2] += area * c_t.y() * c_t.y();
-    covariance[3] += area * c_t.z() * c_t.x();
-    covariance[4] += area * c_t.z() * c_t.y();
-    covariance[5] += area * c_t.z() * c_t.z();
-  }
-
-  // remove sum(area) * (c * transpose(c))
-  covariance[0] -= sum_areas * c.x() * c.x();
-  covariance[1] -= sum_areas * c.y() * c.x();
-  covariance[2] -= sum_areas * c.y() * c.y();
-  covariance[3] -= sum_areas * c.z() * c.x();
-  covariance[4] -= sum_areas * c.z() * c.y();
-  covariance[5] -= sum_areas * c.z() * c.z();
-}
-
-
-// compute the eigen values and vectors of the covariance 
-// matrix and deduces the best linear fitting plane.
-// returns fitting quality
-template < typename K >
-typename K::FT
-fitting_plane_3(const typename K::FT covariance[6], // covariance matrix
-                const typename K::Point_3& c,       // centroid
-                typename K::Plane_3& plane,         // best fit plane
-                const K& )                          // kernel
-{
-  typedef typename K::FT       FT;
-  typedef typename K::Point_3  Point;
-  typedef typename K::Plane_3  Plane;
-  typedef typename K::Vector_3 Vector;
-
-  // solve for eigenvalues and eigenvectors.
-  // eigen values are sorted in descending order, 
-  // eigen vectors are sorted in accordance.
-  FT eigen_values[3];
-  FT eigen_vectors[9];
-  eigen_symmetric<FT>(covariance,3,eigen_vectors,eigen_values);
-
-  // check unicity and build fitting line accordingly
-  if(eigen_values[0] != eigen_values[1] &&
-     eigen_values[0] != eigen_values[2])
-  {
-    // regular case
-    Vector normal(eigen_vectors[6],
-                  eigen_vectors[7],
-                  eigen_vectors[8]);
-    plane = Plane(c,normal);
-    return (FT)1.0 - eigen_values[2] / eigen_values[0];
-  } // end regular case
-  else
-  {
-    // isotropic case (infinite number of directions)
-    // by default: assemble a horizontal plane that goes
-    // through the centroid.
-    plane = Plane(c,Vector(0.0,0.0,1.0));
-    return (FT)0.0;
-  } 
-}
-
-// compute the eigen values and vectors of the covariance 
-// matrix and deduces the best linear fitting line
-// (this is an internal function)
-// returns fitting quality
-template < typename K >
-typename K::FT
-fitting_line_3(const typename K::FT covariance[6], // covariance matrix
-               const typename K::Point_3& c,       // centroid
-               typename K::Line_3& line,           // best fit line
-               const K&  )                         // kernel
-{
-  typedef typename K::FT       FT;
-  typedef typename K::Point_3  Point;
-  typedef typename K::Line_3   Line;
-  typedef typename K::Vector_3 Vector;
-
-  // solve for eigenvalues and eigenvectors.
-  // eigen values are sorted in descending order, 
-  // eigen vectors are sorted in accordance.
-  FT eigen_values[3];
-  FT eigen_vectors[9];
-  eigen_symmetric<FT>(covariance,3,eigen_vectors,eigen_values);
-
-  // check unicity and build fitting line accordingly
-  if(eigen_values[0] != eigen_values[1])
-  {
-    // regular case
-    Vector direction(eigen_vectors[0],eigen_vectors[1],eigen_vectors[2]);
-    line = Line(c,direction);
-    return (FT)1.0 - eigen_values[1] / eigen_values[0];
-  } // end regular case
-  else
-  {
-    // isotropic case (infinite number of directions)
-    // by default: assemble a horizontal plane that goes
-    // through the centroid.
-    line = Line(c,Vector(1.0,0.0,0.0));
-    return (FT)0.0;
-  } 
-}
-
-
-
-// fits a plane to a 3D point 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 plane with default direction)
 template < typename InputIterator, 
-           typename K >
-typename K::FT
-linear_least_squares_fitting_3(InputIterator first,
-                               InputIterator beyond, 
-                               typename K::Plane_3& plane, // best fit plane
-                               typename K::Point_3& c,     // centroid
-                               const K& k,                 // kernel
-                               const typename K::Point_3*) // used for indirection
-{
-  typedef typename K::FT       FT;
-  typedef typename K::Point_3  Point;
-  typedef typename K::Plane_3  Plane;
-  typedef typename K::Vector_3 Vector;
-
-  // precondition: at least one element in the container.
-  CGAL_precondition(first != beyond);
-
-  // compute centroid
-  c = centroid(first,beyond,K());
-
-  // assemble covariance matrix
-  FT covariance[6];
-  assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Point*) NULL);
-
-  // compute fitting plane
-  return fitting_plane_3(covariance,c,plane,k);
-} // end fit plane to point set
-
-// fits a line to a 3D point 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 along x axis)
-template < typename InputIterator, 
-           typename K >
-typename K::FT
-linear_least_squares_fitting_3(InputIterator first,
-                               InputIterator beyond, 
-                               typename K::Line_3& line,  // best fit line
-                               typename K::Point_3& c,    // centroid
-                               const K& k,                // kernel
-                               const typename K::Point_3*)// used for indirection
-{
-  typedef typename K::FT       FT;
-  typedef typename K::Point_3  Point;
-  typedef typename K::Line_3   Line;
-  typedef typename K::Vector_3 Vector;
-
-  // precondition: at least one element in the container.
-  CGAL_precondition(first != beyond);
-
-  // compute centroid
-  c = centroid(first,beyond,K());
-
-  // assemble covariance matrix
-  FT covariance[6];
-  assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Point*) NULL);
-
-  // compute fitting line
-  return fitting_line_3(covariance,c,line,k);
-} // end fit line to point set
-
-// fits a plane to a 3D triangle set
-template < typename InputIterator, 
-           typename K >
-typename K::FT
-linear_least_squares_fitting_3(InputIterator first,
-                               InputIterator beyond, 
-                               typename K::Plane_3& plane,   // best fit plane
-                               typename K::Point_3& c,       // centroid
-                               const K& k,                   // kernel
-                               const typename K::Triangle_3*)// used for indirection
-{
-  typedef typename K::FT          FT;
-  typedef typename K::Point_3     Point;
-  typedef typename K::Vector_3    Vector;
-  typedef typename K::Triangle_3  Triangle;
-
-  // precondition: at least one element in the container.
-  CGAL_precondition(first != beyond);
-
-  // compute centroid
-  c = centroid(first,beyond,K());
-
-  // assemble covariance matrix
-  FT covariance[6];
-  assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Triangle*) NULL);
-
-  // compute fitting plane
-  return fitting_plane_3(covariance,c,plane,k);
-} // end linear_least_squares_fitting_3
-
-// fits a line to a 3D triangle set
-template < typename InputIterator, 
-           typename K >
-typename K::FT
-linear_least_squares_fitting_3(InputIterator first,
-                               InputIterator beyond, 
-                               typename K::Line_3& line,     // best fit line
-                               typename K::Point_3& c,       // centroid
-                               const K& k,                   // kernel
-                               const typename K::Triangle_3*)// used for indirection
-{
-  typedef typename K::FT          FT;
-  typedef typename K::Point_3     Point;
-  typedef typename K::Vector_3    Vector;
-  typedef typename K::Triangle_3  Triangle;
-  typedef typename K::Line_3      Line;
-
-  // precondition: at least one element in the container.
-  CGAL_precondition(first != beyond);
-
-  // compute centroid
-  c = centroid(first,beyond,K());
-
-  // assemble covariance matrix
-  FT covariance[6];
-  assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Triangle*) NULL);
-
-  // compute fitting plane
-  return fitting_line_3(covariance,c,line,k);
-} // end linear_least_squares_fitting_3
-
-} // end namespace CGALi
-
-template < typename InputIterator, 
-           typename K >
+           typename K , typename tag>
 inline
 typename K::FT
 linear_least_squares_fitting_3(InputIterator first,
                                InputIterator beyond, 
                                typename K::Plane_3& plane,
                                typename K::Point_3& centroid,
-                               const K& k)
+                               const K& k,
+			       const tag& t)
 {
   typedef typename std::iterator_traits<InputIterator>::value_type Value_type;
+  //  BOOST_STATIC_ASSERT((boost::is_same<typename CGAL::Algebraic_structure_traits<Value_type>::Algebraic_category,CGAL::Field_with_sqrt_tag>::value));
   return CGALi::linear_least_squares_fitting_3(first, beyond, plane,
-                                               centroid, k, (Value_type*) NULL);
+                                               centroid, k, (Value_type*) NULL, t);
 }
 
 template < typename InputIterator, 
-           typename K >
+           typename K, typename tag >
 inline
 typename K::FT
 linear_least_squares_fitting_3(InputIterator first,
                                InputIterator beyond, 
                                typename K::Line_3& line,
                                typename K::Point_3& centroid,
-                               const K& k)
+                               const K& k,
+			       const tag& t)
 {
   typedef typename std::iterator_traits<InputIterator>::value_type Value_type;
+  //  BOOST_STATIC_ASSERT((boost::is_same<typename CGAL::Algebraic_structure_traits<Value_type>::Algebraic_category,CGAL::Field_with_sqrt_tag>::value));
   return CGALi::linear_least_squares_fitting_3(first, beyond, line,
-                                               centroid, k, (Value_type*) NULL);
+                                               centroid, k, (Value_type*) NULL, t);
 }
 
 template < typename InputIterator, 
-           typename K >
+           typename K, typename tag >
 inline
 typename K::FT
 linear_least_squares_fitting_3(InputIterator first,
                                InputIterator beyond, 
                                typename K::Plane_3& plane,
-                               const K& k)
+                               const K& k,
+			       const tag& t)
 {
   typedef typename std::iterator_traits<InputIterator>::value_type Value_type;
+  //   BOOST_STATIC_ASSERT((boost::is_same<typename CGAL::Algebraic_structure_traits<Value_type>::Algebraic_category,CGAL::Field_with_sqrt_tag>::value));
   typename K::Point_3 centroid;
   return CGALi::linear_least_squares_fitting_3(first, beyond, plane,
-                                               centroid, k,(Value_type*) NULL);
+                                               centroid, k,(Value_type*) NULL, t);
 }
 
 template < typename InputIterator, 
-           typename K >
+           typename K, typename tag >
 inline
 typename K::FT
 linear_least_squares_fitting_3(InputIterator first,
                                InputIterator beyond, 
                                typename K::Line_3& line,
-                               const K& k)
+                               const K& k,
+			       const tag& t)
 {
   typedef typename std::iterator_traits<InputIterator>::value_type Value_type;
+  //   BOOST_STATIC_ASSERT((boost::is_same<typename CGAL::Algebraic_structure_traits<Value_type>::Algebraic_category,CGAL::Field_with_sqrt_tag>::value));
   typename K::Point_3 centroid;
   return CGALi::linear_least_squares_fitting_3(first, beyond, line,
-                                               centroid, k,(Value_type*) NULL);
+                                               centroid, k,(Value_type*) NULL, t);
 }
 
 // deduces the kernel from the points in container.
 template < typename InputIterator, 
            typename Object,
-           typename Point>
+           typename Point, typename tag>
 inline
 typename Kernel_traits<Object>::Kernel::FT
 linear_least_squares_fitting_3(InputIterator first,
                                InputIterator beyond, 
                                Object& object,
-                               Point& centroid)
+                               Point& centroid,
+			       const tag& t)
 {
   typedef typename std::iterator_traits<InputIterator>::value_type Value_type;
+  //   BOOST_STATIC_ASSERT((boost::is_same<typename CGAL::Algebraic_structure_traits<Value_type>::Algebraic_category,CGAL::Field_with_sqrt_tag>::value));
   typedef typename Kernel_traits<Value_type>::Kernel K;
-  return CGAL::linear_least_squares_fitting_3(first,beyond,object,centroid,K());
+  return CGAL::linear_least_squares_fitting_3(first,beyond,object,centroid,K(), t);
 }
 
 // does not return the centroid and deduces the kernel as well.
 template < typename InputIterator, 
-           typename Object >
+           typename Object, typename tag >
 inline
 typename Kernel_traits<Object>::Kernel::FT
 linear_least_squares_fitting_3(InputIterator first,
                                InputIterator beyond, 
-                               Object& object)
+                               Object& object,
+			       const tag& t)
 {
   typedef typename std::iterator_traits<InputIterator>::value_type Value_type;
+  //   BOOST_STATIC_ASSERT((boost::is_same<typename CGAL::Algebraic_structure_traits<Value_type>::Algebraic_category,CGAL::Field_with_sqrt_tag>::value));
   typedef typename Kernel_traits<Value_type>::Kernel K;
-  return CGAL::linear_least_squares_fitting_3(first,beyond,object,K());
+  return CGAL::linear_least_squares_fitting_3(first,beyond,object,K(), t);
 }
 
 CGAL_END_NAMESPACE
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
index 7270e329903..7134d121f48 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
@@ -19,8 +19,6 @@
 #ifndef CGAL_LINEAR_LEAST_SQUARES_FITTING_CIRCLES_2_H
 #define CGAL_LINEAR_LEAST_SQUARES_FITTING_CIRCLES_2_H
 
-#define PI 3.1415926535
-
 #include <CGAL/basic.h>
 #include <CGAL/Object.h>
 #include <CGAL/centroid.h>
@@ -50,7 +48,7 @@ linear_least_squares_fitting_2(InputIterator first,
                                typename K::Point_2& c,     // centroid
                                const K&,                   // kernel
                                const typename K::Circle_2*,// used for indirection
-			       const bool non_standard_geometry)  // true means it is a disk
+			       const CGAL::PCA_dimension_2_tag& tag)
 {
   // types
   typedef typename K::FT       FT;
@@ -63,98 +61,10 @@ linear_least_squares_fitting_2(InputIterator first,
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
-  if(!non_standard_geometry) {
+
+  // ::::::::::DISK:::::::::::::::
   // 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 circles and their combined mass
-  FT mass = 0.0;
-  FT covariance[3] = {0.0,0.0,0.0};
-
-  // assemble 2nd order moment about the origin.  
-  FT temp[4] = {1.0, 0.0,
-		0.0, 1.0};
-  Matrix moment = init_Matrix<K>(2,temp);
-  //  Matrix moment = Matrix(2,true,PI);
-
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    // Now for each circle, construct the 2nd order moment about the origin.
-    // assemble the transformation matrix.
-    const Circle& t = *it;
-
-    // defined for convenience.
-    // FT example = CGAL::to_double(t[0].x());
-    FT radius = std::sqrt(t.squared_radius());
-    FT delta[4] = {radius, 0.0, 
-		   0.0, radius};
-    Matrix transformation = init_Matrix<K>(2,delta);
-    FT area = t.squared_radius();
-    CGAL_assertion(area != 0.0);
-
-    // Find the 2nd order moment for the circle wrt to the origin by an affine transformation.
-    
-    // Transform the standard 2nd order moment using the transformation matrix
-    transformation = area * transformation * moment * LA::transpose(transformation);
-    
-    // Translate the 2nd order moment to the center of the circle.
-    FT x0 = t.center().x();
-    FT y0 = t.center().y();
-
-    // and add to covariance matrix
-    covariance[0] += transformation[0][0] + area * x0*x0;
-    covariance[1] += transformation[0][1] + area * x0*y0;
-    covariance[2] += transformation[1][1] + area * y0*y0;
-
-    mass += area;
-  }
-
-  // Translate the 2nd order moment calculated about the origin to
-  // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.y() * c.y());
-
-  // to remove later
-  std::cout<<covariance[0]<<" "<<covariance[1]<<" "<<covariance[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>(covariance,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;
-  } 
-  }
-  else {
-    // ::::::::::DISK:::::::::::::::
-  // compute centroid
-  c = centroid(first,beyond,K());
+  c = centroid(first,beyond,K(),tag);
 
   // assemble covariance matrix as a semi-definite matrix. 
   // Matrix numbering:
@@ -239,8 +149,117 @@ linear_least_squares_fitting_2(InputIterator first,
     line = Line(c, Vector(1.0, 0.0));
     return (FT)0.0;
   } 
+  
+} // end linear_least_squares_fitting_2 for circle set with 2D tag
+
+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::Circle_2*,// used for indirection
+			        const CGAL::PCA_dimension_1_tag& tag)
+{
+  // 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::Circle_2 Circle;
+  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(),tag);
+
+  // assemble covariance matrix as a semi-definite matrix. 
+  // Matrix numbering:
+  // 0
+  // 1 2
+  //Final combined covariance matrix for all circles and their combined mass
+  FT mass = 0.0;
+  FT covariance[3] = {0.0,0.0,0.0};
+
+  // assemble 2nd order moment about the origin.  
+  FT temp[4] = {1.0, 0.0,
+		0.0, 1.0};
+  Matrix moment = init_Matrix<K>(2,temp);
+
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    // Now for each circle, construct the 2nd order moment about the origin.
+    // assemble the transformation matrix.
+    const Circle& t = *it;
+
+    // defined for convenience.
+    // FT example = CGAL::to_double(t[0].x());
+    FT radius = std::sqrt(t.squared_radius());
+    FT delta[4] = {radius, 0.0, 
+		   0.0, radius};
+    Matrix transformation = init_Matrix<K>(2,delta);
+    FT length = 2 * radius;
+    CGAL_assertion(length != 0.0);
+
+    // Find the 2nd order moment for the circle wrt to the origin by an affine transformation.
+    
+    // Transform the standard 2nd order moment using the transformation matrix
+    transformation = 0.5 * length * transformation * moment * LA::transpose(transformation);
+    
+    // Translate the 2nd order moment to the center of the circle.
+    FT x0 = t.center().x();
+    FT y0 = t.center().y();
+
+    // and add to covariance matrix
+    covariance[0] += transformation[0][0] + length * x0*x0;
+    covariance[1] += transformation[0][1] + length * x0*y0;
+    covariance[2] += transformation[1][1] + length * y0*y0;
+
+    mass += length;
   }
-} // end linear_least_squares_fitting_2 for circle set
+
+  // Translate the 2nd order moment calculated about the origin to
+  // the center of mass to get the covariance.
+  covariance[0] += mass * (-1.0 * c.x() * c.x());
+  covariance[1] += mass * (-1.0 * c.x() * c.y());
+  covariance[2] += mass * (-1.0 * c.y() * c.y());
+
+  // to remove later
+  std::cout<<covariance[0]<<" "<<covariance[1]<<" "<<covariance[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>(covariance,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 circle set with 1D tag
 
 } // end namespace CGALi
 
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_cuboids_3.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_cuboids_3.h
new file mode 100644
index 00000000000..d9456ebd00b
--- /dev/null
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_cuboids_3.h
@@ -0,0 +1,330 @@
+// 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_cuboids_3.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_CUBOIDS_3_H
+#define CGAL_LINEAR_LEAST_SQUARES_FITTING_CUBOIDS_3_H
+
+#include <CGAL/basic.h>
+#include <CGAL/Object.h>
+#include <CGAL/centroid.h>
+#include <CGAL/eigen.h>
+#include <CGAL/util.h>
+
+#include <iterator>
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+// fits a plane to a 3D cuboid set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Plane_3& plane,   // best fit plane
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Iso_cuboid_3*,  // used for indirection
+			       const CGAL::PCA_dimension_3_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Iso_cuboid_3    Iso_cuboid;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  // compute centroid
+  c = centroid(first,beyond,K(),tag);
+
+  // assemble covariance matrix
+  FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
+  assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Iso_cuboid*) NULL,tag);
+
+  // to remove later
+  //  std::cout<<covariance[0]<<" "<<covariance[1]<<" "<<covariance[2]<<" "<<covariance[3]<<" "<<covariance[4]<<" "<<covariance[5]<<std::endl;
+  
+  // compute fitting plane
+  return fitting_plane_3(covariance,c,plane,k);
+
+
+} // end linear_least_squares_fitting_cuboids_3
+
+// fits a plane to a 3D cuboid set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Plane_3& plane,   // best fit plane
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Iso_cuboid_3*,  // used for indirection
+			       const CGAL::PCA_dimension_2_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Iso_cuboid_3    Iso_cuboid;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  // compute centroid
+  c = centroid(first,beyond,K(),tag);
+
+  // assemble covariance matrix
+  FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
+  assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Iso_cuboid*) NULL,tag);
+
+  // to remove later
+  //  std::cout<<covariance[0]<<" "<<covariance[1]<<" "<<covariance[2]<<" "<<covariance[3]<<" "<<covariance[4]<<" "<<covariance[5]<<std::endl;
+  
+  // compute fitting plane
+  return fitting_plane_3(covariance,c,plane,k);
+
+
+} // end linear_least_squares_fitting_cuboids_3
+
+// fits a plane to a 3D cuboid set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Plane_3& plane,   // best fit plane
+                               typename K::Segment_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Iso_cuboid_3*,  // used for indirection
+			       const CGAL::PCA_dimension_1_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Iso_cuboid_3    Iso_cuboid;
+  typedef typename K::Segment_3  Segment;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  std::list<Segment> segments;
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Iso_cuboid& t = *it;
+    segments.push_back(t[0],t[1]);
+    segments.push_back(t[1],t[2]);
+    segments.push_back(t[2],t[3]);    
+    segments.push_back(t[3],t[0]);
+    segments.push_back(t[4],t[5]);
+    segments.push_back(t[5],t[6]);
+    segments.push_back(t[6],t[7]);
+    segments.push_back(t[7],t[4]);
+    segments.push_back(t[5],t[0]);
+    segments.push_back(t[1],t[6]);
+    segments.push_back(t[2],t[7]);
+    segments.push_back(t[3],t[4]);
+  }
+
+  // compute fitting plane
+  return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,k,(Segment*)NULL,tag);
+
+} // end linear_least_squares_fitting_cuboids_3
+
+// fits a plane to a 3D cuboid set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Plane_3& plane,   // best fit plane
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Iso_cuboid_3*,  // used for indirection
+			       const CGAL::PCA_dimension_0_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Iso_cuboid_3    Iso_cuboid;
+  typedef typename K::Point_3  Point;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  std::list<Point> points;
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Iso_cuboid& t = *it;
+    points.push_back(t[0]);
+    points.push_back(t[1]);
+    points.push_back(t[2]);    
+    points.push_back(t[3]);
+  }
+
+  // compute fitting plane
+  return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,k,(Point*)NULL,tag);
+
+} // end linear_least_squares_fitting_cuboids_3
+
+// fits a line to a 3D cuboid set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Line_3& line,     // best fit line
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Iso_cuboid_3*,  // used for indirection
+			       const CGAL::PCA_dimension_3_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Iso_cuboid_3    Iso_cuboid;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  // compute centroid
+  c = centroid(first,beyond,K(),tag);
+
+  // assemble covariance matrix
+  FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
+  assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Iso_cuboid*) NULL,tag);
+  
+  // compute fitting line
+  return fitting_line_3(covariance,c,line,k);
+  
+} // end linear_least_squares_fitting_cuboids_3
+
+// fits a line to a 3D cuboid set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Line_3& line,   // best fit line
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Iso_cuboid_3*,  // used for indirection
+			       const CGAL::PCA_dimension_2_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Iso_cuboid_3    Iso_cuboid;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  // compute centroid
+  c = centroid(first,beyond,K(),tag);
+
+  // assemble covariance matrix
+  FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
+  assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Iso_cuboid*) NULL,tag);
+
+  // to remove later
+  //  std::cout<<covariance[0]<<" "<<covariance[1]<<" "<<covariance[2]<<" "<<covariance[3]<<" "<<covariance[4]<<" "<<covariance[5]<<std::endl;
+  
+  // compute fitting line
+  return fitting_line_3(covariance,c,line,k);
+
+
+} // end linear_least_squares_fitting_cuboids_3
+
+// fits a line to a 3D cuboid set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Line_3& line,   // best fit line
+                               typename K::Segment_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Iso_cuboid_3*,  // used for indirection
+			       const CGAL::PCA_dimension_1_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Iso_cuboid_3    Iso_cuboid;
+  typedef typename K::Segment_3  Segment;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  std::list<Segment> segments;
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Iso_cuboid& t = *it;
+    segments.push_back(t[0],t[1]);
+    segments.push_back(t[1],t[2]);
+    segments.push_back(t[2],t[3]);    
+    segments.push_back(t[3],t[0]);
+    segments.push_back(t[4],t[5]);
+    segments.push_back(t[5],t[6]);
+    segments.push_back(t[6],t[7]);
+    segments.push_back(t[7],t[4]);
+    segments.push_back(t[5],t[0]);
+    segments.push_back(t[1],t[6]);
+    segments.push_back(t[2],t[7]);
+    segments.push_back(t[3],t[4]);
+  }
+
+  // compute fitting line
+  return linear_least_squares_fitting_3(segments.begin(),segments.end(),line,c,k,(Segment*)NULL,tag);
+
+} // end linear_least_squares_fitting_cuboids_3
+
+// fits a line to a 3D cuboid set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Line_3& line,   // best fit line
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Iso_cuboid_3*,  // used for indirection
+			       const CGAL::PCA_dimension_0_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Iso_cuboid_3    Iso_cuboid;
+  typedef typename K::Point_3  Point;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  std::list<Point> points;
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Iso_cuboid& t = *it;
+    points.push_back(t[0]);
+    points.push_back(t[1]);
+    points.push_back(t[2]);    
+    points.push_back(t[3]);
+  }
+
+  // compute fitting line
+  return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,k,(Point*)NULL,tag);
+
+} // end linear_least_squares_fitting_cuboids_3
+
+} // end namespace CGALi
+
+CGAL_END_NAMESPACE
+
+#endif // CGAL_LINEAR_LEAST_SQUARES_FITTING_CUBOIDS_3_H
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_points_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_points_2.h
index a0d16ab95c5..20810a10bda 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_points_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_points_2.h
@@ -45,7 +45,7 @@ linear_least_squares_fitting_2(InputIterator first,
                                typename K::Point_2& c,     // centroid
                                const K&,                   // kernel
                                const typename K::Point_2*,// used for indirection
-			       const bool non_standard_geometry)  // not useful 
+			       const CGAL::PCA_dimension_0_tag& t)   
 {
   // types
   typedef typename K::FT       FT;
@@ -59,7 +59,7 @@ linear_least_squares_fitting_2(InputIterator first,
   CGAL_precondition(first != beyond);
 
   // compute centroid
-  c = centroid(first,beyond,K());
+  c = centroid(first,beyond,K(),t);
 
   // assemble covariance matrix as a semi-definite matrix. 
   // Matrix numbering:
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_points_3.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_points_3.h
new file mode 100644
index 00000000000..f50f21c2434
--- /dev/null
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_points_3.h
@@ -0,0 +1,105 @@
+// 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_points_3.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_POINTS_3_H
+#define CGAL_LINEAR_LEAST_SQUARES_FITTING_POINTS_3_H
+
+#include <CGAL/basic.h>
+#include <CGAL/Object.h>
+#include <CGAL/centroid.h>
+#include <CGAL/eigen.h>
+#include <CGAL/util.h>
+
+#include <iterator>
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+// fits a plane to a 3D point 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 plane with default direction)
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Plane_3& plane, // best fit plane
+                               typename K::Point_3& c,     // centroid
+                               const K& k,                 // kernel
+                               const typename K::Point_3*, // used for indirection
+			       const CGAL::PCA_dimension_0_tag& tag)
+{
+  typedef typename K::FT       FT;
+  typedef typename K::Point_3  Point;
+  typedef typename K::Plane_3  Plane;
+  typedef typename K::Vector_3 Vector;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+
+  // compute centroid
+  c = centroid(first,beyond,K(),tag);
+
+  // assemble covariance matrix
+  FT covariance[6];
+  assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Point*) NULL,tag);
+
+  // compute fitting plane
+  return fitting_plane_3(covariance,c,plane,k);
+} // end fit plane to point set
+
+// fits a line to a 3D point 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 along x axis)
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Line_3& line,  // best fit line
+                               typename K::Point_3& c,    // centroid
+                               const K& k,                // kernel
+                               const typename K::Point_3*, // used for indirection
+			       const CGAL::PCA_dimension_0_tag& tag)
+{
+  typedef typename K::FT       FT;
+  typedef typename K::Point_3  Point;
+  typedef typename K::Line_3   Line;
+  typedef typename K::Vector_3 Vector;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+
+  // compute centroid
+  c = centroid(first,beyond,K(),tag);
+
+  // assemble covariance matrix
+  FT covariance[6];
+  assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Point*) NULL,tag);
+
+  // compute fitting line
+  return fitting_line_3(covariance,c,line,k);
+} // end fit line to point set
+
+} // end namespace CGALi
+
+CGAL_END_NAMESPACE
+
+#endif // CGAL_LINEAR_LEAST_SQUARES_FITTING_POINTS_3_H
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
index 0fc7d38541a..106c24838f6 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
@@ -48,7 +48,7 @@ linear_least_squares_fitting_2(InputIterator first,
                                typename K::Point_2& c,     // centroid
                                const K&,                   // kernel
                                const typename K::Iso_rectangle_2*,// used for indirection
-			       const bool non_standard_geometry)  // true means it is an empty rectangle
+			       const CGAL::PCA_dimension_2_tag& tag)
 {
   // types
   typedef typename K::FT       FT;
@@ -62,9 +62,9 @@ linear_least_squares_fitting_2(InputIterator first,
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
-  if(!non_standard_geometry) {
+
   // compute centroid
-  c = centroid(first,beyond,K());
+  c = centroid(first,beyond,K(),tag);
 
   // assemble covariance matrix as a semi-definite matrix. 
   // Matrix numbering:
@@ -125,7 +125,7 @@ linear_least_squares_fitting_2(InputIterator first,
   covariance[2] += mass * (-1.0 * c.y() * c.y());
 
   // to remove later
-  std::cout<<covariance[0]<<" "<<covariance[1]<<" "<<covariance[2]<<std::endl;
+  //  std::cout<<covariance[0]<<" "<<covariance[1]<<" "<<covariance[2]<<std::endl;
 
   // solve for eigenvalues and eigenvectors.
   // eigen values are sorted in descending order, 
@@ -153,24 +153,75 @@ linear_least_squares_fitting_2(InputIterator first,
     line = Line(c, Vector(1.0, 0.0));
     return (FT)0.0;
   } 
-  }
-  else {
-    std::list<Segment_2> segments;
-    
-    for(InputIterator it = first;
-	it != beyond;
-	it++)
-    {
-      const Iso_rectangle& t = *it;
-      segments.push_back(Segment_2(t[0],t[1]));
-      segments.push_back(Segment_2(t[1],t[2]));
-      segments.push_back(Segment_2(t[2],t[3]));      
-      segments.push_back(Segment_2(t[3],t[0]));      
-    }    
+} // end linear_least_squares_fitting_2 for rectangle set with 2D tag
 
-    return linear_least_squares_fitting_2(segments.begin(),segments.end(),line,c,K());
-  }
-} // end linear_least_squares_fitting_2 for rectangle set
+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::Iso_rectangle_2*,// used for indirection
+			       const CGAL::PCA_dimension_1_tag& tag)
+{
+  // types
+  typedef typename K::Iso_rectangle_2 Iso_rectangle;
+  typedef typename K::Segment_2         Segment_2;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+
+  std::list<Segment_2> segments;
+  
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Iso_rectangle& t = *it;
+    segments.push_back(Segment_2(t[0],t[1]));
+    segments.push_back(Segment_2(t[1],t[2]));
+    segments.push_back(Segment_2(t[2],t[3]));      
+    segments.push_back(Segment_2(t[3],t[0]));      
+  }    
+
+  return linear_least_squares_fitting_2(segments.begin(),segments.end(),line,c,K(),tag);
+
+} // end linear_least_squares_fitting_2 for rectangle set with 1D tag
+
+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::Iso_rectangle_2*,// used for indirection
+			       const CGAL::PCA_dimension_0_tag& tag)
+{
+  // types
+  typedef typename K::Iso_rectangle_2 Iso_rectangle;
+  typedef typename K::Point_2         Point_2;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+
+  std::list<Point_2> points;
+  
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Iso_rectangle& t = *it;
+    points.push_back(Point_2(t[0]));
+    points.push_back(Point_2(t[1]));
+    points.push_back(Point_2(t[2]));      
+    points.push_back(Point_2(t[3]));      
+  }    
+
+  return linear_least_squares_fitting_2(points.begin(),points.end(),line,c,K(),tag);
+
+} // end linear_least_squares_fitting_2 for rectangle set with 0D tag
 
 } // end namespace CGALi
 
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
index 8c3f3ef5df0..6ca55dbf163 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
@@ -48,7 +48,7 @@ linear_least_squares_fitting_2(InputIterator first,
                                typename K::Point_2& c,     // centroid
                                const K&,                   // kernel
                                const typename K::Segment_2*,// used for indirection
-			       const bool non_standard_geometry)  // not useful 
+			       const CGAL::PCA_dimension_1_tag& tag)   
 {
   // types
   typedef typename K::FT       FT;
@@ -63,7 +63,7 @@ linear_least_squares_fitting_2(InputIterator first,
   CGAL_precondition(first != beyond);
 
   // compute centroid
-  c = centroid(first,beyond,K());
+  c = centroid(first,beyond,K(),tag);
   // assemble covariance matrix as a semi-definite matrix. 
   // Matrix numbering:
   // 0
@@ -111,6 +111,7 @@ linear_least_squares_fitting_2(InputIterator first,
   covariance[0] += mass * (-1.0 * c.x() * c.x());
   covariance[1] += mass * (-1.0 * c.x() * c.y());
   covariance[2] += mass * (-1.0 * c.y() * c.y());
+  std::cout<<covariance[0]<<" "<<covariance[1]<<" "<<covariance[2]<<" "<<(std::sqrt(2)/3.0)<<std::endl;
 
   // to remove later
   //  std::cout<<covariance[0]<<" "<<covariance[1]<<" "<<covariance[2]<<" "<<(std::sqrt(2)/3.0)<<std::endl;
@@ -141,7 +142,37 @@ linear_least_squares_fitting_2(InputIterator first,
     line = Line(c, Vector(1.0, 0.0));
     return (FT)0.0;
   } 
-} // end linear_least_squares_fitting_2 for segment set
+} // end linear_least_squares_fitting_2 for segment set with 1D tag
+
+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& k,                   // kernel
+                               const typename K::Segment_2*,// used for indirection
+			       const CGAL::PCA_dimension_0_tag& tag)   
+{
+  // types
+  typedef typename K::Point_2  Point;
+  typedef typename K::Segment_2 Segment;
+ 
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+
+  std::list<Point> points;  
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Segment& s = *it;
+    points.push_back(s[0]);
+    points.push_back(s[1]);
+  } 
+  return linear_least_squares_fitting_2(points.begin(),points.end(),line,c,k,(Point*)NULL,tag);
+
+} // end linear_least_squares_fitting_2 for segment set with 1D tag
 
 } // end namespace CGALi
 
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_3.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_3.h
new file mode 100644
index 00000000000..9dd790d2913
--- /dev/null
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_3.h
@@ -0,0 +1,172 @@
+// 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_segments_3.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_3_H
+#define CGAL_LINEAR_LEAST_SQUARES_FITTING_SEGMENTS_3_H
+
+#include <CGAL/basic.h>
+#include <CGAL/Object.h>
+#include <CGAL/centroid.h>
+#include <CGAL/eigen.h>
+#include <CGAL/util.h>
+
+#include <iterator>
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+// fits a plane to a 3D segment set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Plane_3& plane,   // best fit plane
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Segment_3*,  // used for indirection
+			       const CGAL::PCA_dimension_1_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Segment_3  Segment;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  // compute centroid
+  c = centroid(first,beyond,K(),tag);
+
+  // assemble covariance matrix
+  FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
+  assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Segment*) NULL,tag);
+
+  // to remove later
+  std::cout<<covariance[0]<<" "<<covariance[1]<<" "<<covariance[2]<<" "<<covariance[3]<<" "<<covariance[4]<<" "<<covariance[5]<<std::endl;
+  
+  // compute fitting plane
+  return fitting_plane_3(covariance,c,plane,k);
+
+} // end linear_least_squares_fitting_segments_3
+
+// fits a plane to a 3D segment set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Plane_3& plane,   // best fit plane
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Segment_3*,  // used for indirection
+			       const CGAL::PCA_dimension_0_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Segment_3  Segment;
+  typedef typename K::Point_3  Point;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  std::list<Point> points;
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Segment& t = *it;
+    points.push_back(t[0]);
+    points.push_back(t[1]);
+  }
+
+  // compute fitting plane
+  return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,k,(Point*)NULL,tag);
+
+} // end linear_least_squares_fitting_segments_3
+
+// fits a line to a 3D segment set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Line_3& line,     // best fit line
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Segment_3*,  // used for indirection
+			       const CGAL::PCA_dimension_1_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Segment_3  Segment;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  // compute centroid
+  c = centroid(first,beyond,K(),tag);
+
+  // assemble covariance matrix
+  FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
+  assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Segment*) NULL,tag);
+
+  // to remove later
+  std::cout<<covariance[0]<<" "<<covariance[1]<<" "<<covariance[2]<<" "<<covariance[3]<<" "<<covariance[4]<<" "<<covariance[5]<<std::endl;
+
+  // compute fitting line
+  return fitting_line_3(covariance,c,line,k);
+  
+} // end linear_least_squares_fitting_segments_3
+
+// fits a plane to a 3D segment set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Line_3& line,   // best fit plane
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Segment_3*,  // used for indirection
+			       const CGAL::PCA_dimension_0_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Segment_3  Segment;
+  typedef typename K::Point_3  Point;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  std::list<Point> points;
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Segment& t = *it;
+    points.push_back(t[0]);
+    points.push_back(t[1]);
+  }
+
+  // compute fitting plane
+  return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,k,(Point*)NULL,tag);
+
+} // end linear_least_squares_fitting_segments_3
+
+} // end namespace CGALi
+
+CGAL_END_NAMESPACE
+
+#endif // CGAL_LINEAR_LEAST_SQUARES_FITTING_SEGMENTS_3_H
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_spheres_3.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_spheres_3.h
new file mode 100644
index 00000000000..55c87413b83
--- /dev/null
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_spheres_3.h
@@ -0,0 +1,164 @@
+// 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_spheres_3.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_SPHERES_3_H
+#define CGAL_LINEAR_LEAST_SQUARES_FITTING_SPHERES_3_H
+
+#include <CGAL/basic.h>
+#include <CGAL/Object.h>
+#include <CGAL/centroid.h>
+#include <CGAL/eigen.h>
+#include <CGAL/util.h>
+
+#include <iterator>
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+// fits a plane to a 3D sphere set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Plane_3& plane,   // best fit plane
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Sphere_3*,  // used for indirection
+			       const CGAL::PCA_dimension_3_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Sphere_3    Sphere;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  // compute centroid
+  c = centroid(first,beyond,K(),tag);
+
+  // assemble covariance matrix
+  FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
+  assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Sphere*) NULL,tag);
+  
+  // compute fitting plane
+  return fitting_plane_3(covariance,c,plane,k);
+
+} // end linear_least_squares_fitting_spheres_3
+
+// fits a plane to a 3D sphere set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Plane_3& plane,   // best fit plane
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Sphere_3*,  // used for indirection
+			       const CGAL::PCA_dimension_2_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Sphere_3    Sphere;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  // compute centroid
+  c = centroid(first,beyond,K(),tag);
+
+  // assemble covariance matrix
+  FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
+  assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Sphere*) NULL,tag);
+  
+  // compute fitting plane
+  return fitting_plane_3(covariance,c,plane,k);
+
+} // end linear_least_squares_fitting_spheres_3
+
+// fits a line to a 3D sphere set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Line_3& line,     // best fit line
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Sphere_3*,  // used for indirection
+			       const CGAL::PCA_dimension_3_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Sphere_3  Sphere;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  // compute centroid
+  c = centroid(first,beyond,K(),tag);
+  
+  // assemble covariance matrix
+  FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
+  assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Sphere*) NULL,tag);
+
+  // to remove later
+  //  std::cout<<covariance[0]<<" "<<covariance[1]<<" "<<covariance[2]<<" "<<covariance[3]<<" "<<covariance[4]<<" "<<covariance[5]<<std::endl;
+
+  // compute fitting line
+  return fitting_line_3(covariance,c,line,k);
+  
+} // end linear_least_squares_fitting_spheres_3
+
+// fits a line to a 3D sphere set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Line_3& line,     // best fit line
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Sphere_3*,  // used for indirection
+			       const CGAL::PCA_dimension_2_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Sphere_3  Sphere;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+
+  // compute centroid
+  c = centroid(first,beyond,K(),tag);
+  
+  // assemble covariance matrix
+  FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
+  assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Sphere*) NULL,tag);
+
+  // to remove later
+  //  std::cout<<covariance[0]<<" "<<covariance[1]<<" "<<covariance[2]<<" "<<covariance[3]<<" "<<covariance[4]<<" "<<covariance[5]<<std::endl;
+
+  // compute fitting line
+  return fitting_line_3(covariance,c,line,k);
+
+} // end linear_least_squares_fitting_spheres_3
+
+} // end namespace CGALi
+
+CGAL_END_NAMESPACE
+
+#endif // CGAL_LINEAR_LEAST_SQUARES_FITTING_SPHERES_3_H
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_tetrahedrons_3.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_tetrahedrons_3.h
new file mode 100644
index 00000000000..8062d9639b6
--- /dev/null
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_tetrahedrons_3.h
@@ -0,0 +1,325 @@
+// Copyright (c) 2005  INRIA Sophia-Antipolis (France).
+// All rights reserved.
+//
+// This file is part of CGAL (woidww.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_tetrahedrons_3.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_TETRAHEDRONS_3_H
+#define CGAL_LINEAR_LEAST_SQUARES_FITTING_TETRAHEDRONS_3_H
+
+#include <CGAL/basic.h>
+#include <CGAL/Object.h>
+#include <CGAL/centroid.h>
+#include <CGAL/eigen.h>
+#include <CGAL/util.h>
+
+#include <iterator>
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+// fits a plane to a 3D tetrahedron set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Plane_3& plane,   // best fit plane
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Tetrahedron_3*,  // used for indirection
+			       const CGAL::PCA_dimension_3_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Tetrahedron_3    Tetrahedron;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  // compute centroid
+  c = centroid(first,beyond,K(),tag);
+
+  // assemble covariance matrix
+  FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
+  assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Tetrahedron*) NULL,tag);
+
+  // to remove later
+  std::cout<<covariance[0]<<" "<<covariance[1]<<" "<<covariance[2]<<" "<<covariance[3]<<" "<<covariance[4]<<" "<<covariance[5]<<std::endl;
+  
+  // compute fitting plane
+  return fitting_plane_3(covariance,c,plane,k);
+
+
+} // end linear_least_squares_fitting_tetrahedrons_3
+
+// fits a plane to a 3D tetrahedron set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Plane_3& plane,   // best fit plane
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Tetrahedron_3*,  // used for indirection
+			       const CGAL::PCA_dimension_2_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Tetrahedron_3    Tetrahedron;
+  typedef typename K::Triangle_3 Triangle;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+
+  std::list<Triangle> triangles;
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Tetrahedron& t = *it;
+    triangles.push_back(t[0],t[1],t[2]);
+    triangles.push_back(t[0],t[2],t[3]);
+    triangles.push_back(t[0],t[3],t[1]); 
+    triangles.push_back(t[3],t[1],t[2]);
+ }
+
+  // compute fitting plane
+  return linear_least_squares_fitting_3(triangles.begin(),triangles.end(),plane,c,k,(Triangle*)NULL,tag);
+
+} // end linear_least_squares_fitting_tetrahedrons_3
+
+// fits a plane to a 3D tetrahedron set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Plane_3& plane,   // best fit plane
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Tetrahedron_3*,  // used for indirection
+			       const CGAL::PCA_dimension_1_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Tetrahedron_3    Tetrahedron;
+  typedef typename K::Segment_3 Segment;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+
+  std::list<Segment> segments;
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Tetrahedron& t = *it;
+    segments.push_back(t[0],t[1]);
+    segments.push_back(t[1],t[2]);
+    segments.push_back(t[1],t[3]);
+    segments.push_back(t[2],t[3]); 
+    segments.push_back(t[0],t[2]);
+    segments.push_back(t[0],t[3]);
+ }
+
+  // compute fitting plane
+  return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,k,(Segment*)NULL,tag);
+
+} // end linear_least_squares_fitting_tetrahedrons_3
+
+// fits a plane to a 3D tetrahedron set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Plane_3& plane,   // best fit plane
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Tetrahedron_3*,  // used for indirection
+			       const CGAL::PCA_dimension_0_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Tetrahedron_3    Tetrahedron;
+  typedef typename K::Point_3 Point;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+
+  std::list<Point> points;
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Tetrahedron& t = *it;
+    points.push_back(t[0]);
+    points.push_back(t[1]);
+    points.push_back(t[2]); 
+    points.push_back(t[3]);
+ }
+
+  // compute fitting plane
+  return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,k,(Point*)NULL,tag);
+
+} // end linear_least_squares_fitting_tetrahedrons_3
+
+// fits a line to a 3D tetrahedron set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Line_3& line,     // best fit line
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Tetrahedron_3*,  // used for indirection
+			       const CGAL::PCA_dimension_3_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Tetrahedron_3    Tetrahedron;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  // compute centroid
+  c = centroid(first,beyond,K(),tag);
+
+  // assemble covariance matrix
+  FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
+  assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Tetrahedron*) NULL,tag);
+  
+  // to remove later
+  std::cout<<covariance[0]<<" "<<covariance[1]<<" "<<covariance[2]<<" "<<covariance[3]<<" "<<covariance[4]<<" "<<covariance[5]<<std::endl;
+  
+  // compute fitting line
+  return fitting_line_3(covariance,c,line,k);
+  
+} // end linear_least_squares_fitting_tetrahedrons_3
+
+// fits a line to a 3D tetrahedron set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Line_3& line,   // best fit line
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Tetrahedron_3*,  // used for indirection
+			       const CGAL::PCA_dimension_2_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Tetrahedron_3    Tetrahedron;
+  typedef typename K::Triangle_3 Triangle;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+
+  std::list<Triangle> triangles;
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Tetrahedron& t = *it;
+    triangles.push_back(t[0],t[1],t[2]);
+    triangles.push_back(t[0],t[2],t[3]);
+    triangles.push_back(t[0],t[3],t[1]); 
+    triangles.push_back(t[3],t[1],t[2]);
+ }
+
+  // compute fitting line
+  return linear_least_squares_fitting_3(triangles.begin(),triangles.end(),line,c,k,(Triangle*)NULL,tag);
+
+} // end linear_least_squares_fitting_tetrahedrons_3
+
+// fits a line to a 3D tetrahedron set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Line_3& line,   // best fit line
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Tetrahedron_3*,  // used for indirection
+			       const CGAL::PCA_dimension_1_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Tetrahedron_3    Tetrahedron;
+  typedef typename K::Segment_3 Segment;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+
+  std::list<Segment> segments;
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Tetrahedron& t = *it;
+    segments.push_back(t[0],t[1]);
+    segments.push_back(t[1],t[2]);
+    segments.push_back(t[1],t[3]);
+    segments.push_back(t[2],t[3]); 
+    segments.push_back(t[0],t[2]);
+    segments.push_back(t[0],t[3]);
+ }
+
+  // compute fitting line
+  return linear_least_squares_fitting_3(segments.begin(),segments.end(),line,c,k,(Segment*)NULL,tag);
+
+} // end linear_least_squares_fitting_tetrahedrons_3
+
+// fits a line to a 3D tetrahedron set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Line_3& line,   // best fit line
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Tetrahedron_3*,  // used for indirection
+			       const CGAL::PCA_dimension_0_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Tetrahedron_3    Tetrahedron;
+  typedef typename K::Point_3 Point;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+
+  std::list<Point> points;
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Tetrahedron& t = *it;
+    points.push_back(t[0]);
+    points.push_back(t[1]);
+    points.push_back(t[2]); 
+    points.push_back(t[3]);
+ }
+
+  // compute fitting line
+  return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,k,(Point*)NULL,tag);
+
+} // end linear_least_squares_fitting_tetrahedrons_3
+
+} // end namespace CGALi
+
+CGAL_END_NAMESPACE
+
+#endif // CGAL_LINEAR_LEAST_SQUARES_FITTING_TETRAHEDRONS_3_H
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
index 6c30e1c2cc1..68fda7f6dc7 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
@@ -49,7 +49,7 @@ linear_least_squares_fitting_2(InputIterator first,
                                typename K::Point_2& c,     // centroid
                                const K&,                   // kernel
                                const typename K::Triangle_2*,// used for indirection
-			       const bool non_standard_geometry)  // true means it is a hollow triangle
+			       const CGAL::PCA_dimension_2_tag& tag)
 {
   // types
   typedef typename K::FT       FT;
@@ -64,9 +64,8 @@ linear_least_squares_fitting_2(InputIterator first,
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  if(!non_standard_geometry) {
   // compute centroid
-  c = centroid(first,beyond,K());
+  c = centroid(first,beyond,K(),tag);
 
   // assemble covariance matrix as a semi-definite matrix. 
   // Matrix numbering:
@@ -97,24 +96,24 @@ linear_least_squares_fitting_2(InputIterator first,
 		   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);
+    FT area = 0.5 * std::abs(LA::determinant(transformation));
+    CGAL_assertion(area!=0);
 
     // Find the 2nd order moment for the triangle wrt to the origin by an affine transformation.
     
     // Transform the standard 2nd order moment using the transformation matrix
-    transformation = detA * transformation * moment * LA::transpose(transformation);
+    transformation = 2 * area * transformation * moment * LA::transpose(transformation);
     
     // Translate the 2nd order moment to (x0,y0).
     FT xav0 = (delta[0]+delta[1])/3.0;
     FT yav0 = (delta[2]+delta[3])/3.0;
     
     // and add to the covariance matrix
-    covariance[0] += transformation[0][0] + 0.5 * detA * (x0*xav0*2 + x0*x0);
-    covariance[1] += transformation[0][1] + 0.5 * detA * (x0*yav0 + xav0*y0 + x0*y0);
-    covariance[2] += transformation[1][1] + 0.5 * detA * (y0*yav0*2 + y0*y0);
+    covariance[0] += transformation[0][0] + area * (x0*xav0*2 + x0*x0);
+    covariance[1] += transformation[0][1] + area * (x0*yav0 + xav0*y0 + x0*y0);
+    covariance[2] += transformation[1][1] + area * (y0*yav0*2 + y0*y0);
 
-    mass += 0.5 *detA;
+    mass += area;
   }
   
   // Translate the 2nd order moment calculated about the origin to 
@@ -122,6 +121,8 @@ linear_least_squares_fitting_2(InputIterator first,
   covariance[0] += mass * (-1.0 * c.x() * c.x());
   covariance[1] += mass * (-1.0 * c.x() * c.y());
   covariance[2] += mass * (-1.0 * c.y() * c.y());
+  // to remove later
+  std::cout<<covariance[0]<<" "<<covariance[1]<<" "<<covariance[2]<<std::endl;
 
   // solve for eigenvalues and eigenvectors.
   // eigen values are sorted in descending order, 
@@ -149,23 +150,73 @@ linear_least_squares_fitting_2(InputIterator first,
     line = Line(c, Vector(1.0, 0.0));
     return (FT)0.0;
   } 
-  }
-  else {
-    std::list<Segment_2> segments;
-    
-    for(InputIterator it = first;
-	it != beyond;
-	it++)
-    {
-      const Triangle& t = *it;
-      segments.push_back(Segment_2(t[0],t[1]));
-      segments.push_back(Segment_2(t[1],t[2]));
-      segments.push_back(Segment_2(t[2],t[0]));      
-    }    
+} // end linear_least_squares_fitting_2 for triangle set with 2D tag
 
-    return linear_least_squares_fitting_2(segments.begin(),segments.end(),line,c,K());
-  }
-} // end linear_least_squares_fitting_2 for triangle set
+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
+			       const CGAL::PCA_dimension_1_tag& tag)
+{
+  // types
+
+  typedef typename K::Triangle_2 Triangle;
+  typedef typename K::Segment_2  Segment_2;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  std::list<Segment_2> segments;  
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Triangle& t = *it;
+    segments.push_back(Segment_2(t[0],t[1]));
+    segments.push_back(Segment_2(t[1],t[2]));
+    segments.push_back(Segment_2(t[2],t[0]));      
+  }    
+  
+  return linear_least_squares_fitting_2(segments.begin(),segments.end(),line,c,K(),tag);
+  
+} // end linear_least_squares_fitting_2 for triangle set with 1D tag
+
+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
+			       const CGAL::PCA_dimension_0_tag& tag)
+{
+  // types
+
+  typedef typename K::Triangle_2 Triangle;
+  typedef typename K::Point_2  Point_2;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  std::list<Point_2> points;  
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Triangle& t = *it;
+    points.push_back(Point_2(t[0]));
+    points.push_back(Point_2(t[1]));
+    points.push_back(Point_2(t[2]));      
+  }    
+  
+  return linear_least_squares_fitting_2(points.begin(),points.end(),line,c,K(),tag);
+  
+} // end linear_least_squares_fitting_2 for triangle set with 0D tag
 
 } // end namespace CGALi
 
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h
new file mode 100644
index 00000000000..f2cb1d6d347
--- /dev/null
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h
@@ -0,0 +1,244 @@
+// 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_3.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_3_H
+#define CGAL_LINEAR_LEAST_SQUARES_FITTING_TRIANGLES_3_H
+
+#include <CGAL/basic.h>
+#include <CGAL/Object.h>
+#include <CGAL/centroid.h>
+#include <CGAL/eigen.h>
+#include <CGAL/util.h>
+
+#include <iterator>
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+// fits a plane to a 3D triangle set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Plane_3& plane,   // best fit plane
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Triangle_3*,  // used for indirection
+			       const CGAL::PCA_dimension_2_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Triangle_3  Triangle;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  // compute centroid
+  c = centroid(first,beyond,K(),tag);
+
+  // assemble covariance matrix
+  FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
+  assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Triangle*) NULL,tag);
+
+  // to remove later
+  std::cout<<covariance[0]<<" "<<covariance[1]<<" "<<covariance[2]<<" "<<covariance[3]<<" "<<covariance[4]<<" "<<covariance[5]<<std::endl;
+  
+  // compute fitting plane
+  return fitting_plane_3(covariance,c,plane,k);
+
+} // end linear_least_squares_fitting_triangles_3
+
+// fits a plane to a 3D triangle set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Plane_3& plane,   // best fit plane
+                               typename K::Segment_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Triangle_3*,  // used for indirection
+			       const CGAL::PCA_dimension_1_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Triangle_3  Triangle;
+  typedef typename K::Segment_3  Segment;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  std::list<Segment> segments;
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Triangle& t = *it;
+    segments.push_back(t[0],t[1]);
+    segments.push_back(t[1],t[2]);
+    segments.push_back(t[2],t[0]);    
+  }
+
+  // compute fitting plane
+  return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,k,(Segment*)NULL,tag);
+
+} // end linear_least_squares_fitting_triangles_3
+
+// fits a plane to a 3D triangle set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Plane_3& plane,   // best fit plane
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Triangle_3*,  // used for indirection
+			       const CGAL::PCA_dimension_0_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Triangle_3  Triangle;
+  typedef typename K::Point_3  Point;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  std::list<Point> points;
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Triangle& t = *it;
+    points.push_back(t[0]);
+    points.push_back(t[1]);
+    points.push_back(t[2]);    
+  }
+
+  // compute fitting plane
+  return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,k,(Point*)NULL,tag);
+
+} // end linear_least_squares_fitting_triangles_3
+
+// fits a line to a 3D triangle set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Line_3& line,     // best fit line
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Triangle_3*,  // used for indirection
+			       const CGAL::PCA_dimension_2_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Triangle_3  Triangle;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  // compute centroid
+  c = centroid(first,beyond,K(),tag);
+
+  // assemble covariance matrix
+  FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
+  assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Triangle*) NULL,tag);
+
+  // to remove later
+  std::cout<<covariance[0]<<" "<<covariance[1]<<" "<<covariance[2]<<" "<<covariance[3]<<" "<<covariance[4]<<" "<<covariance[5]<<std::endl;
+
+  // compute fitting line
+  return fitting_line_3(covariance,c,line,k);
+  
+} // end linear_least_squares_fitting_triangles_3
+
+// fits a line to a 3D triangle set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Line_3& line,   // best fit line
+                               typename K::Segment_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Triangle_3*,  // used for indirection
+			       const CGAL::PCA_dimension_1_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Triangle_3  Triangle;
+  typedef typename K::Segment_3  Segment;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  std::list<Segment> segments;
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Triangle& t = *it;
+    segments.push_back(t[0],t[1]);
+    segments.push_back(t[1],t[2]);
+    segments.push_back(t[2],t[0]);    
+  }
+
+  // compute fitting line
+  return linear_least_squares_fitting_3(segments.begin(),segments.end(),line,c,k,(Segment*)NULL,tag);
+
+} // end linear_least_squares_fitting_triangles_3
+
+// fits a line to a 3D triangle set
+template < typename InputIterator, 
+           typename K >
+typename K::FT
+linear_least_squares_fitting_3(InputIterator first,
+                               InputIterator beyond, 
+                               typename K::Line_3& line,   // best fit line
+                               typename K::Point_3& c,       // centroid
+                               const K& k,                   // kernel
+                               const typename K::Triangle_3*,  // used for indirection
+			       const CGAL::PCA_dimension_0_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Triangle_3  Triangle;
+  typedef typename K::Point_3  Point;
+
+  // precondition: at least one element in the container.
+  CGAL_precondition(first != beyond);
+  
+  std::list<Point> points;
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Triangle& t = *it;
+    points.push_back(t[0]);
+    points.push_back(t[1]);
+    points.push_back(t[2]);    
+  }
+
+  // compute fitting line
+  return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,k,(Point*)NULL,tag);
+
+} // end linear_least_squares_fitting_triangles_3
+
+} // end namespace CGALi
+
+CGAL_END_NAMESPACE
+
+#endif // CGAL_LINEAR_LEAST_SQUARES_FITTING_TRIANGLES_3_H
diff --git a/Principal_component_analysis/include/CGAL/util.h b/Principal_component_analysis/include/CGAL/util.h
index 04da2b95674..ddb65843ea2 100644
--- a/Principal_component_analysis/include/CGAL/util.h
+++ b/Principal_component_analysis/include/CGAL/util.h
@@ -22,6 +22,7 @@
 #include <CGAL/basic.h>
 #include <CGAL/Object.h>
 #include <CGAL/Linear_algebraCd.h>
+#include <CGAL/PCA_tags.h>
 
 CGAL_BEGIN_NAMESPACE
 
@@ -42,6 +43,685 @@ init_Matrix(int n, typename K::FT entries[]) {
   return ret;
 } // end initialization of matrix
 
+// assemble covariance matrix from a point set 
+template < typename InputIterator,
+           typename K >
+void
+assemble_covariance_matrix_3(InputIterator first,
+                             InputIterator beyond, 
+                             typename K::FT covariance[6], // covariance matrix
+                             const typename K::Point_3& c, // centroid
+                             const K& ,                    // kernel
+                             const typename K::Point_3*,   // used for indirection
+			     const CGAL::PCA_dimension_0_tag& tag)
+{
+  typedef typename K::FT       FT;
+  typedef typename K::Point_3  Point;
+  typedef typename K::Vector_3 Vector;
+
+  // Matrix numbering:
+  // 0          
+  // 1 2
+  // 3 4 5          
+  covariance[0] = covariance[1] = covariance[2] = 
+  covariance[3] = covariance[4] = covariance[5] = (FT)0.0;
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Point& p = *it;
+    Vector d = p - c;
+    covariance[0] += d.x() * d.x();
+    covariance[1] += d.x() * d.y();
+    covariance[2] += d.y() * d.y();
+    covariance[3] += d.x() * d.z();
+    covariance[4] += d.y() * d.z();
+    covariance[5] += d.z() * d.z();
+  }
+}
+
+// assemble covariance matrix from a triangle set 
+template < typename InputIterator,
+           typename K >
+void
+assemble_covariance_matrix_3(InputIterator first,
+                             InputIterator beyond, 
+                             typename K::FT covariance[6], // covariance matrix
+                             const typename K::Point_3& c, // centroid
+                             const K& ,                    // kernel
+                             const typename K::Triangle_3*,// used for indirection
+			     const CGAL::PCA_dimension_2_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Point_3     Point;
+  typedef typename K::Vector_3    Vector;
+  typedef typename K::Triangle_3  Triangle;
+  typedef typename CGAL::Linear_algebraCd<FT> LA;
+  typedef typename LA::Matrix Matrix;
+
+  // assemble covariance matrix as a semi-definite matrix. 
+  // Matrix numbering:
+  // 0
+  // 1 2 
+  // 3 4 5
+  //Final combined covariance matrix for all triangles and their combined mass
+  FT mass = 0.0;
+
+  // assemble 2nd order moment about the origin.  
+  FT temp[9] = {1.0/12.0, 1.0/24.0, 1.0/24.0,
+		1.0/24.0, 1.0/12.0, 1.0/24.0,
+                1.0/24.0, 1.0/24.0, 1.0/12.0};
+  Matrix moment = init_Matrix<K>(3,temp);
+
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    // Now for each triangle, construct the 2nd order moment about the origin.
+    // assemble the transformation matrix.
+    const Triangle& t = *it;
+
+    // defined for convenience.
+    // FT example = CGAL::to_double(t[0].x());
+    FT delta[9] = {t[0].x(), t[1].x(), t[2].x(), 
+		   t[0].y(), t[1].y(), t[2].y(),
+                   t[0].z(), t[1].z(), t[2].z()};
+    Matrix transformation = init_Matrix<K>(3,delta);
+    FT area = std::sqrt(t.squared_area());
+    CGAL_assertion(area != 0.0);
+
+    // Find the 2nd order moment for the triangle wrt to the origin by an affine transformation.
+    
+    // Transform the standard 2nd order moment using the transformation matrix
+    transformation = 2 * area * transformation * moment * LA::transpose(transformation);
+    
+    // and add to covariance matrix
+    covariance[0] += transformation[0][0];
+    covariance[1] += transformation[1][0];
+    covariance[2] += transformation[1][1];
+    covariance[3] += transformation[2][0];
+    covariance[4] += transformation[2][1];
+    covariance[5] += transformation[2][2];
+
+    mass += area;
+  }
+
+  // Translate the 2nd order moment calculated about the origin to
+  // the center of mass to get the covariance.
+  covariance[0] += mass * (-1.0 * c.x() * c.x());
+  covariance[1] += mass * (-1.0 * c.x() * c.y());
+  covariance[2] += mass * (-1.0 * c.y() * c.y());
+  covariance[3] += mass * (-1.0 * c.z() * c.x());
+  covariance[4] += mass * (-1.0 * c.z() * c.y());
+  covariance[5] += mass * (-1.0 * c.z() * c.z());
+
+}
+
+// assemble covariance matrix from a cuboid set 
+template < typename InputIterator,
+           typename K >
+void
+assemble_covariance_matrix_3(InputIterator first,
+                             InputIterator beyond, 
+                             typename K::FT covariance[6], // covariance matrix
+                             const typename K::Point_3& c, // centroid
+                             const K& ,                    // kernel
+                             const typename K::Iso_cuboid_3*,// used for indirection
+			     const CGAL::PCA_dimension_3_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Point_3     Point;
+  typedef typename K::Vector_3    Vector;
+  typedef typename K::Iso_cuboid_3    Iso_cuboid;
+  typedef typename CGAL::Linear_algebraCd<FT> LA;
+  typedef typename LA::Matrix Matrix;
+
+  // assemble covariance matrix as a semi-definite matrix. 
+  // Matrix numbering:
+  // 0
+  // 1 2 
+  // 3 4 5
+  //Final combined covariance matrix for all cuboids and their combined mass
+  FT mass = 0.0;
+
+  // assemble 2nd order moment about the origin.  
+  FT temp[9] = {1.0/3.0, 1.0/4.0, 1.0/4.0,
+		1.0/4.0, 1.0/3.0, 1.0/4.0,
+                1.0/4.0, 1.0/4.0, 1.0/3.0};
+  Matrix moment = init_Matrix<K>(3,temp);
+
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    // Now for each cuboid, construct the 2nd order moment about the origin.
+    // assemble the transformation matrix.
+    const Iso_cuboid& t = *it;
+
+    // defined for convenience.
+    // FT example = CGAL::to_double(t[0].x());
+    FT x0 = t[0].x();
+    FT y0 = t[0].y();
+    FT z0 = t[0].z();
+    FT delta[9] = {t[1].x()-x0, t[3].x()-x0, t[5].x()-x0, 
+		   t[1].y()-y0, t[3].y()-y0, t[5].y()-y0,
+                   t[1].z()-z0, t[3].z()-z0, t[5].z()-z0};
+    Matrix transformation = init_Matrix<K>(3,delta);
+    FT volume = t.volume();
+    CGAL_assertion(volume != 0.0);
+
+    // Find the 2nd order moment for the cuboid wrt to the origin by an affine transformation.
+    
+    // Transform the standard 2nd order moment using the transformation matrix
+    transformation = volume * transformation * moment * LA::transpose(transformation);
+    
+    // Translate the 2nd order moment to the minimum corner (x0,y0,z0) of the cuboid.
+    FT xav0 = (delta[0] + delta[1] + delta[2])/4.0;
+    FT yav0 = (delta[3] + delta[4] + delta[5])/4.0;
+    FT zav0 = (delta[6] + delta[7] + delta[8])/4.0;
+
+    // and add to covariance matrix
+    covariance[0] += transformation[0][0] + volume * (2*x0*xav0 + x0*x0);
+    covariance[1] += transformation[1][0] + volume * (xav0*y0 + yav0*x0 + x0*y0);
+    covariance[2] += transformation[1][1] + volume * (2*y0*yav0 + y0*y0);
+    covariance[3] += transformation[2][0] + volume * (x0*zav0 + xav0*z0 + x0*z0);
+    covariance[4] += transformation[2][1] + volume * (yav0*z0 + y0*zav0 + z0*y0);
+    covariance[5] += transformation[2][2] + volume * (2*zav0*z0 + z0*z0);
+
+    mass += volume;
+  }
+
+  // Translate the 2nd order moment calculated about the origin to
+  // the center of mass to get the covariance.
+  covariance[0] += mass * (-1.0 * c.x() * c.x());
+  covariance[1] += mass * (-1.0 * c.x() * c.y());
+  covariance[2] += mass * (-1.0 * c.y() * c.y());
+  covariance[3] += mass * (-1.0 * c.z() * c.x());
+  covariance[4] += mass * (-1.0 * c.z() * c.y());
+  covariance[5] += mass * (-1.0 * c.z() * c.z());
+
+}
+
+// assemble covariance matrix from a cuboid set 
+template < typename InputIterator,
+           typename K >
+void
+assemble_covariance_matrix_3(InputIterator first,
+                             InputIterator beyond, 
+                             typename K::FT covariance[6], // covariance matrix
+                             const typename K::Point_3& c, // centroid
+                             const K& ,                    // kernel
+                             const typename K::Iso_cuboid_3*,// used for indirection
+			     const CGAL::PCA_dimension_2_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Point_3     Point;
+  typedef typename K::Vector_3    Vector;
+  typedef typename K::Iso_cuboid_3    Iso_cuboid;
+  typedef typename CGAL::Linear_algebraCd<FT> LA;
+  typedef typename LA::Matrix Matrix;
+
+  // assemble covariance matrix as a semi-definite matrix. 
+  // Matrix numbering:
+  // 0
+  // 1 2 
+  // 3 4 5
+  //Final combined covariance matrix for all cuboids and their combined mass
+  FT mass = 0.0;
+
+  // assemble 2nd order moment about the origin.  
+  FT temp[9] = {7.0/3.0, 1.5, 1.5,
+		1.5, 7.0/3.0, 1.5,
+                1.5, 1.5, 7.0/3.0};
+  Matrix moment = init_Matrix<K>(3,temp);
+
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    // Now for each cuboid, construct the 2nd order moment about the origin.
+    // assemble the transformation matrix.
+    const Iso_cuboid& t = *it;
+
+    // defined for convenience.
+    // FT example = CGAL::to_double(t[0].x());
+    FT x0 = t[0].x();
+    FT y0 = t[0].y();
+    FT z0 = t[0].z();
+    FT delta[9] = {t[1].x()-x0, t[3].x()-x0, t[5].x()-x0, 
+		   t[1].y()-y0, t[3].y()-y0, t[5].y()-y0,
+                   t[1].z()-z0, t[3].z()-z0, t[5].z()-z0};
+    Matrix transformation = init_Matrix<K>(3,delta);
+    FT volume = t.volume();
+    CGAL_assertion(volume != 0.0);
+
+    // Find the 2nd order moment for the cuboid wrt to the origin by an affine transformation.
+    
+    // Transform the standard 2nd order moment using the transformation matrix
+    transformation = volume * transformation * moment * LA::transpose(transformation);
+    
+    // Translate the 2nd order moment to the minimum corner (x0,y0,z0) of the cuboid.
+    FT xav0 = (delta[0] + delta[1] + delta[2])/4.0;
+    FT yav0 = (delta[3] + delta[4] + delta[5])/4.0;
+    FT zav0 = (delta[6] + delta[7] + delta[8])/4.0;
+
+    // and add to covariance matrix
+    covariance[0] += transformation[0][0] + volume * (2*x0*xav0 + x0*x0);
+    covariance[1] += transformation[1][0] + volume * (xav0*y0 + yav0*x0 + x0*y0);
+    covariance[2] += transformation[1][1] + volume * (2*y0*yav0 + y0*y0);
+    covariance[3] += transformation[2][0] + volume * (x0*zav0 + xav0*z0 + x0*z0);
+    covariance[4] += transformation[2][1] + volume * (yav0*z0 + y0*zav0 + z0*y0);
+    covariance[5] += transformation[2][2] + volume * (2*zav0*z0 + z0*z0);
+
+    mass += volume;
+  }
+
+  // Translate the 2nd order moment calculated about the origin to
+  // the center of mass to get the covariance.
+  covariance[0] += mass * (-1.0 * c.x() * c.x());
+  covariance[1] += mass * (-1.0 * c.x() * c.y());
+  covariance[2] += mass * (-1.0 * c.y() * c.y());
+  covariance[3] += mass * (-1.0 * c.z() * c.x());
+  covariance[4] += mass * (-1.0 * c.z() * c.y());
+  covariance[5] += mass * (-1.0 * c.z() * c.z());
+
+}
+
+// assemble covariance matrix from a sphere set 
+template < typename InputIterator,
+           typename K >
+void
+assemble_covariance_matrix_3(InputIterator first,
+                             InputIterator beyond, 
+                             typename K::FT covariance[6], // covariance matrix
+                             const typename K::Point_3& c, // centroid
+                             const K& ,                    // kernel
+                             const typename K::Sphere_3*,// used for indirection
+			     const CGAL::PCA_dimension_3_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Point_3     Point;
+  typedef typename K::Vector_3    Vector;
+  typedef typename K::Sphere_3  Sphere;
+  typedef typename CGAL::Linear_algebraCd<FT> LA;
+  typedef typename LA::Matrix Matrix;
+
+  // assemble covariance matrix as a semi-definite matrix. 
+  // Matrix numbering:
+  // 0
+  // 1 2 
+  // 3 4 5
+  //Final combined covariance matrix for all spheres and their combined mass
+  FT mass = 0.0;
+
+  // assemble 2nd order moment about the origin.  
+  FT temp[9] = {4.0/15.0, 0.0,     0.0,
+		0.0,     4.0/15.0, 0.0,
+                0.0,     0.0,     4.0/15.0};
+  Matrix moment = init_Matrix<K>(3,temp);
+
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    // Now for each sphere, construct the 2nd order moment about the origin.
+    // assemble the transformation matrix.
+    const Sphere& t = *it;
+
+    // defined for convenience.
+    // FT example = CGAL::to_double(t[0].x());
+    FT radius = std::sqrt(t.squared_radius());
+    FT delta[9] = {radius, 0.0, 0.0, 
+		   0.0, radius, 0.0,
+                   0.0, 0.0, radius};
+    Matrix transformation = init_Matrix<K>(3,delta);
+    FT area = 4/3.0 * radius*t.squared_radius();
+    CGAL_assertion(area != 0.0);
+
+    // Find the 2nd order moment for the sphere wrt to the origin by an affine transformation.
+    
+    // Transform the standard 2nd order moment using the transformation matrix
+    transformation = (3.0/4.0) * area * transformation * moment * LA::transpose(transformation);
+    
+    // Translate the 2nd order moment to the center of the sphere.
+    FT x0 = t.center().x();
+    FT y0 = t.center().y();
+    FT z0 = t.center().z();
+
+    // and add to covariance matrix
+    covariance[0] += transformation[0][0] + area * x0*x0;
+    covariance[1] += transformation[1][0] + area * x0*y0;
+    covariance[2] += transformation[1][1] + area * y0*y0;
+    covariance[3] += transformation[2][0] + area * x0*z0;
+    covariance[4] += transformation[2][1] + area * z0*y0;
+    covariance[5] += transformation[2][2] + area * z0*z0;
+
+    mass += area;
+  }
+
+  // Translate the 2nd order moment calculated about the origin to
+  // the center of mass to get the covariance.
+  covariance[0] += mass * (-1.0 * c.x() * c.x());
+  covariance[1] += mass * (-1.0 * c.x() * c.y());
+  covariance[2] += mass * (-1.0 * c.y() * c.y());
+  covariance[3] += mass * (-1.0 * c.z() * c.x());
+  covariance[4] += mass * (-1.0 * c.z() * c.y());
+  covariance[5] += mass * (-1.0 * c.z() * c.z());
+
+}
+
+// assemble covariance matrix from a tetrahedron set 
+template < typename InputIterator,
+           typename K >
+void
+assemble_covariance_matrix_3(InputIterator first,
+                             InputIterator beyond, 
+                             typename K::FT covariance[6], // covariance matrix
+                             const typename K::Point_3& c, // centroid
+                             const K& ,                    // kernel
+                             const typename K::Tetrahedron_3*,// used for indirection
+			     const CGAL::PCA_dimension_3_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Point_3     Point;
+  typedef typename K::Vector_3    Vector;
+  typedef typename K::Tetrahedron_3  Tetrahedron;
+  typedef typename CGAL::Linear_algebraCd<FT> LA;
+  typedef typename LA::Matrix Matrix;
+
+  // assemble covariance matrix as a semi-definite matrix. 
+  // Matrix numbering:
+  // 0
+  // 1 2 
+  // 3 4 5
+  //Final combined covariance matrix for all tetrahedrons and their combined mass
+  FT mass = 0.0;
+
+  // assemble 2nd order moment about the origin.  
+  FT temp[9] = {1.0/60.0,  1.0/120.0, 1.0/120.0,
+		1.0/120.0, 1.0/60.0,  1.0/120.0,
+		1.0/120.0, 1.0/120.0, 1.0/60.0};
+  Matrix moment = init_Matrix<K>(3,temp);
+
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    // Now for each tetrahedron, construct the 2nd order moment about the origin.
+    // assemble the transformation matrix.
+    const Tetrahedron& t = *it;
+
+    // defined for convenience.
+    // FT example = CGAL::to_double(t[0].x());
+    FT x0 = t[0].x();
+    FT y0 = t[0].y();
+    FT z0 = t[0].z();
+
+    FT delta[9] = {t[1].x()-x0, t[2].x()-x0, t[3].x()-x0, 
+		   t[1].y()-y0, t[2].y()-y0, t[3].y()-y0,
+                   t[1].z()-z0, t[2].z()-z0, t[3].z()-z0};
+    Matrix transformation = init_Matrix<K>(3,delta);
+    FT volume = t.volume();
+    CGAL_assertion(volume != 0.0);
+
+    // Find the 2nd order moment for the tetrahedron wrt to the origin by an affine transformation.
+    
+    // Transform the standard 2nd order moment using the transformation matrix
+    transformation = 6 * volume * transformation * moment * LA::transpose(transformation);
+    
+    // Translate the 2nd order moment to the center of the tetrahedron.
+    FT xav0 = (delta[0]+delta[1]+delta[2])/4.0;
+    FT yav0 = (delta[3]+delta[4]+delta[5])/4.0;
+    FT zav0 = (delta[6]+delta[7]+delta[8])/4.0;
+
+    // and add to covariance matrix
+    covariance[0] += transformation[0][0] + volume * (2*x0*xav0 + x0*x0);
+    covariance[1] += transformation[1][0] + volume * (xav0*y0 + yav0*x0 + x0*y0);
+    covariance[2] += transformation[1][1] + volume * (2*y0*yav0 + y0*y0);
+    covariance[3] += transformation[2][0] + volume * (x0*zav0 + xav0*z0 + x0*z0);
+    covariance[4] += transformation[2][1] + volume * (yav0*z0 + y0*zav0 + z0*y0);
+    covariance[5] += transformation[2][2] + volume * (2*zav0*z0 + z0*z0);
+
+    mass += volume;
+  }
+
+  // Translate the 2nd order moment calculated about the origin to
+  // the center of mass to get the covariance.
+  covariance[0] += mass * (-1.0 * c.x() * c.x());
+  covariance[1] += mass * (-1.0 * c.x() * c.y());
+  covariance[2] += mass * (-1.0 * c.y() * c.y());
+  covariance[3] += mass * (-1.0 * c.z() * c.x());
+  covariance[4] += mass * (-1.0 * c.z() * c.y());
+  covariance[5] += mass * (-1.0 * c.z() * c.z());
+
+}
+
+// assemble covariance matrix from a segment set 
+template < typename InputIterator,
+           typename K >
+void
+assemble_covariance_matrix_3(InputIterator first,
+                             InputIterator beyond, 
+                             typename K::FT covariance[6], // covariance matrix
+                             const typename K::Point_3& c, // centroid
+                             const K& ,                    // kernel
+                             const typename K::Segment_3*,// used for indirection
+			     const CGAL::PCA_dimension_1_tag& tag)
+{
+  typedef typename K::FT          FT;
+  typedef typename K::Point_3     Point;
+  typedef typename K::Vector_3    Vector;
+  typedef typename K::Segment_3  Segment;
+  typedef typename CGAL::Linear_algebraCd<FT> LA;
+  typedef typename LA::Matrix Matrix;
+
+  // assemble covariance matrix as a semi-definite matrix. 
+  // Matrix numbering:
+  // 0
+  // 1 2 
+  // 3 4 5
+  //Final combined covariance matrix for all segments and their combined mass
+  FT mass = 0.0;
+
+  // assemble 2nd order moment about the origin.  
+  FT temp[9] = {1.0, 0.5, 0.0,
+		0.5, 1.0, 0.0,
+                0.0, 0.0, 0.0};
+  Matrix moment = 1.0/3.0 * init_Matrix<K>(3,temp);
+
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    // Now for each segment, construct the 2nd order moment about the origin.
+    // assemble the transformation matrix.
+    const Segment& t = *it;
+
+    // defined for convenience.
+    // FT example = CGAL::to_double(t[0].x());
+    FT delta[9] = {t[0].x(), t[1].x(), 0.0, 
+		   t[0].y(), t[1].y(), 0.0,
+                   t[0].z(), t[1].z(), 1.0};
+    Matrix transformation = init_Matrix<K>(3,delta);
+    FT length = std::sqrt(t.squared_length());
+    CGAL_assertion(length != 0.0);
+
+    // Find the 2nd order moment for the segment wrt to the origin by an affine transformation.
+    
+    // Transform the standard 2nd order moment using the transformation matrix
+    transformation = length * transformation * moment * LA::transpose(transformation);
+  std::cout<<covariance[0]<<" "<<covariance[1]<<" "<<covariance[2]<<" "<<covariance[3]<<" "<<covariance[4]<<" "<<covariance[5]<<std::endl;    
+    // and add to covariance matrix
+    covariance[0] += transformation[0][0];
+    covariance[1] += transformation[1][0];
+    covariance[2] += transformation[1][1];
+    covariance[3] += transformation[2][0];
+    covariance[4] += transformation[2][1];
+    covariance[5] += transformation[2][2];
+  std::cout<<covariance[0]<<" "<<covariance[1]<<" "<<covariance[2]<<" "<<covariance[3]<<" "<<covariance[4]<<" "<<covariance[5]<<std::endl;
+    mass += length;
+  }
+
+  // Translate the 2nd order moment calculated about the origin to
+  // the center of mass to get the covariance.
+  covariance[0] += mass * (-1.0 * c.x() * c.x());
+  covariance[1] += mass * (-1.0 * c.x() * c.y());
+  covariance[2] += mass * (-1.0 * c.y() * c.y());
+  covariance[3] += mass * (-1.0 * c.z() * c.x());
+  covariance[4] += mass * (-1.0 * c.z() * c.y());
+  covariance[5] += mass * (-1.0 * c.z() * c.z());
+
+}
+
+/*
+// assemble covariance matrix from a triangle set 
+template < typename InputIterator,
+           typename K >
+void
+assemble_covariance_matrix_3_pierre(InputIterator first,
+                             InputIterator beyond, 
+                             typename K::FT covariance[6], // covariance matrix
+                             const typename K::Point_3& c, // centroid
+                             const K& ,                    // kernel
+                             const typename K::Triangle_3*)// used for indirection
+{
+  typedef typename K::FT       FT;
+  typedef typename K::Point_3  Point;
+  typedef typename K::Vector_3 Vector;
+  typedef typename K::Triangle_3 Triangle;
+
+  // Matrix numbering:
+  // 0          
+  // 1 2
+  // 3 4 5          
+  covariance[0] = covariance[1] = covariance[2] = 
+  covariance[3] = covariance[4] = covariance[5] = (FT)0.0;
+  FT sum_areas = 0.0;
+  for(InputIterator it = first;
+      it != beyond;
+      it++)
+  {
+    const Triangle& triangle = *it;
+    FT area = std::sqrt(triangle.squared_area());
+    Point c_t = K().construct_centroid_3_object()(triangle[0],triangle[1],triangle[2]);
+    sum_areas += area;
+
+    // e1 = ab, e2 = ac
+    Vector e1 = Vector(triangle[0],triangle[1]);
+    Vector e2 = Vector(triangle[0],triangle[2]);
+
+    FT c1 = (FT)(2.0 * area * 10.0 / 72.0);
+    FT c2 = (FT)(2.0 * area *  7.0 / 72.0);
+        
+    covariance[0] += c1*(e1[0]*e1[0] + e2[0]*e2[0]) + (FT)2.0*c2*e1[0]*e2[0];
+    covariance[1] += c1*(e1[1]*e1[0] + e2[1]*e2[0]) + c2*(e1[1]*e2[0] + e1[0]*e2[1]);
+    covariance[2] += c1*(e1[1]*e1[1] + e2[1]*e2[1]) + (FT)2.0*c2*e1[1]*e2[1];
+    covariance[3] += c1*(e1[2]*e1[0] + e2[2]*e2[0]) + c2*(e1[2]*e2[0] + e1[0]*e2[2]);
+    covariance[4] += c1*(e1[2]*e1[1] + e2[2]*e2[1]) + c2*(e1[2]*e2[1] + e1[1]*e2[2]);
+    covariance[5] += c1*(e1[2]*e1[2] + e2[2]*e2[2]) + (FT)2.0*c2*e1[2]*e2[2];
+    
+    // add area(t) c(t) * transpose(c(t))
+    covariance[0] += area * c_t.x() * c_t.x();
+    covariance[1] += area * c_t.y() * c_t.x();
+    covariance[2] += area * c_t.y() * c_t.y();
+    covariance[3] += area * c_t.z() * c_t.x();
+    covariance[4] += area * c_t.z() * c_t.y();
+    covariance[5] += area * c_t.z() * c_t.z();
+  }
+
+  // remove sum(area) * (c * transpose(c))
+  covariance[0] -= sum_areas * c.x() * c.x();
+  covariance[1] -= sum_areas * c.y() * c.x();
+  covariance[2] -= sum_areas * c.y() * c.y();
+  covariance[3] -= sum_areas * c.z() * c.x();
+  covariance[4] -= sum_areas * c.z() * c.y();
+  covariance[5] -= sum_areas * c.z() * c.z();
+}
+*/
+// compute the eigen values and vectors of the covariance 
+// matrix and deduces the best linear fitting plane.
+// returns fitting quality
+template < typename K >
+typename K::FT
+fitting_plane_3(const typename K::FT covariance[6], // covariance matrix
+                const typename K::Point_3& c,       // centroid
+                typename K::Plane_3& plane,         // best fit plane
+                const K& )                          // kernel
+{
+  typedef typename K::FT       FT;
+  typedef typename K::Point_3  Point;
+  typedef typename K::Plane_3  Plane;
+  typedef typename K::Vector_3 Vector;
+
+  // solve for eigenvalues and eigenvectors.
+  // eigen values are sorted in descending order, 
+  // eigen vectors are sorted in accordance.
+  FT eigen_values[3];
+  FT eigen_vectors[9];
+  eigen_symmetric<FT>(covariance,3,eigen_vectors,eigen_values);
+
+  // check unicity and build fitting line accordingly
+  if(eigen_values[0] != eigen_values[1] &&
+     eigen_values[0] != eigen_values[2])
+  {
+    // regular case
+    Vector normal(eigen_vectors[6],
+                  eigen_vectors[7],
+                  eigen_vectors[8]);
+    plane = Plane(c,normal);
+    return (FT)1.0 - eigen_values[2] / eigen_values[0];
+  } // end regular case
+  else
+  {
+    // isotropic case (infinite number of directions)
+    // by default: assemble a horizontal plane that goes
+    // through the centroid.
+    plane = Plane(c,Vector(0.0,0.0,1.0));
+    return (FT)0.0;
+  } 
+}
+
+// compute the eigen values and vectors of the covariance 
+// matrix and deduces the best linear fitting line
+// (this is an internal function)
+// returns fitting quality
+template < typename K >
+typename K::FT
+fitting_line_3(const typename K::FT covariance[6], // covariance matrix
+               const typename K::Point_3& c,       // centroid
+               typename K::Line_3& line,           // best fit line
+               const K&  )                         // kernel
+{
+  typedef typename K::FT       FT;
+  typedef typename K::Point_3  Point;
+  typedef typename K::Line_3   Line;
+  typedef typename K::Vector_3 Vector;
+
+  // solve for eigenvalues and eigenvectors.
+  // eigen values are sorted in descending order, 
+  // eigen vectors are sorted in accordance.
+  FT eigen_values[3];
+  FT eigen_vectors[9];
+  eigen_symmetric<FT>(covariance,3,eigen_vectors,eigen_values);
+
+  // check unicity and build fitting line accordingly
+  if(eigen_values[0] != eigen_values[1])
+  {
+    // regular case
+    Vector direction(eigen_vectors[0],eigen_vectors[1],eigen_vectors[2]);
+    line = Line(c,direction);
+    return (FT)1.0 - eigen_values[1] / eigen_values[0];
+  } // end regular case
+  else
+  {
+    // isotropic case (infinite number of directions)
+    // by default: assemble a horizontal plane that goes
+    // through the centroid.
+    line = Line(c,Vector(1.0,0.0,0.0));
+    return (FT)0.0;
+  } 
+}
+
 } // end namespace CGALi
 
 CGAL_END_NAMESPACE