turn Vec3List into a class and give it some convenience methods

modules/geom/pymod/export_vec3.cc

@@ -18,8 +18,9 @@
 //------------------------------------------------------------------------------
 #include <boost/python.hpp>
 #include <boost/python/suite/indexing/vector_indexing_suite.hpp>
-
+#include <ost/geom/vec3.hh>
 #include <ost/geom/geom.hh>
+#include <ost/export_helper/vector.hh>
 
 using namespace boost::python;
 
@@ -63,10 +64,15 @@ void export_Vec3()
     .add_property("y", &Vec3::GetY, &Vec3::SetY)
     .add_property("z", &Vec3::GetZ, &Vec3::SetZ)
   ;
+  
   def("Normalize", &NormalizeV3);
   def("Cross", &Cross);
   
   class_<Vec3List>("Vec3List", init<>())
     .def(vector_indexing_suite<Vec3List>())
+    .def(ost::VectorAdditions<Vec3List>())
+    .add_property("center", &Vec3List::GetCenter)
+    .add_property("inertia", &Vec3List::GetInertia)
+    .add_property("principal_axes", &Vec3List::GetPrincipalAxes)
   ;
 }

modules/geom/src/CMakeLists.txt

@@ -24,6 +24,7 @@ composite2_op.hh
 composite3_op.hh
 aligned_cuboid.hh
 quat.hh
+
 point_cloud.hh
 module_config.hh
 )
@@ -32,6 +33,7 @@ set(OST_GEOM_SOURCES
 mat2.cc
 mat3.cc
 mat4.cc
+vec3.cc
 vecmat2_op.cc
 vecmat3_op.cc
 vecmat4_op.cc

modules/geom/src/vec3.cc

+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2010 by the OpenStructure authors
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation; either version 3.0 of the License, or (at your option)
+// any later version.
+// This library is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with this library; if not, write to the Free Software Foundation, Inc.,
+// 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+//------------------------------------------------------------------------------
+#include <algorithm>
+
+#include <Eigen/SVD>
+#include "vec3.hh"
+
+
+namespace geom {
+
+typedef Eigen::Matrix3f EMat3;
+
+Mat3 Vec3List::GetInertia() const
+{
+  Mat3 cov(0,0,0,0,0,0,0,0,0);
+  Vec3 center=this->GetCenter();
+  for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+    Vec3 p=*i-center;
+    cov(0,0)+=p.y*p.y+p.z*p.z;
+    cov(1,1)+=p.x*p.x+p.z*p.z;
+    cov(2,2)+=p.x*p.x+p.y*p.y;
+    cov(0,1)-=p.x*p.y;
+    cov(1,2)-=p.y*p.z;
+    cov(0,2)-=p.x*p.z;
+  }  
+  cov(1,0)=cov(0,1);    
+  cov(2,1)=cov(1,2);    
+  cov(2,0)=cov(0,2);  
+  return cov;
+}
+
+Mat3 Vec3List::GetPrincipalAxes() const
+{
+  Mat3 inertia=this->GetInertia();  
+  EMat3 inertia_mat(*reinterpret_cast<EMat3*>(&inertia));
+
+  Eigen::SVD<EMat3> svd(inertia_mat);
+  EMat3 rot=svd.matrixU();
+  Mat3 axes(*reinterpret_cast<Mat3*>(&rot));
+  return axes;
+}
+
+Vec3 Vec3List::GetCenter() const
+{
+  Vec3 center;
+  if (this->empty()) {
+    return center;
+  }
+  for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+    center+=*i;
+  }
+  return center/=this->size();
+}
+
+
+}

modules/geom/src/vec3.hh

@@ -192,10 +192,32 @@ inline std::ostream& operator<<(std::ostream& os, const Vec3& v)
 
 #include <ost/geom/vec2.hh>
 #include <ost/geom/vec4.hh>
+#include <ost/geom/mat3.hh>
 
 namespace geom {
 
-typedef std::vector<Vec3> Vec3List;
+class Vec3List : public std::vector<Vec3> {
+public:
+  typedef std::vector<Vec3> base_type;
+  Vec3List() : base_type() {}
+  
+  Vec3List(size_t size, const Vec3& value=Vec3()) : base_type(size, value) {}
+  Vec3List(base_type::iterator b, base_type::iterator e): base_type(b, e) { }
+  
+  Vec3List(const Vec3List& rhs) : base_type(rhs) { }
+  
+  Vec3List& operator=(const Vec3List& rhs) 
+  {
+    *this=rhs;
+    return *this;
+  }
+  Mat3 GetInertia() const;
+  
+  Vec3 GetCenter() const;
+  
+  Mat3 GetPrincipalAxes() const;
+};
+
 
 inline Vec3::Vec3(const Vec2& v): x(v.x), y(v.y), z(0.0) { }
 inline Vec3::Vec3(const Vec4& v): x(v.x/v.w), y(v.y/v.w), z(v.z/v.w) { }