[Cctbx-cvs] scitbx/include/scitbx/math eigensystem.h, 1.6, 1.7 principal_axes_of_inertia.h, 1.8, 1.9

[Nicholas K. S. <nks...@us...>]

Update of /cvsroot/cctbx/scitbx/include/scitbx/math
In directory sc8-pr-cvs4.sourceforge.net:/tmp/cvs-serv11881/include/scitbx/math

Modified Files:
	eigensystem.h principal_axes_of_inertia.h 
Log Message:
Infrastructure for 2D linear algebra.  Added classes scitbx::vec2,
scitbx::mat2, and scitbx::sym_mat2.  In the eigensystem namespace, a
new constructor real_symmetric(scitbx::sym_mat2) is added.  A new
class, principal_axes_of_inertia_2d now handles the 2D case.


Index: eigensystem.h
===================================================================
RCS file: /cvsroot/cctbx/scitbx/include/scitbx/math/eigensystem.h,v
retrieving revision 1.6
retrieving revision 1.7
diff -C2 -d -r1.6 -r1.7
*** eigensystem.h	21 Jan 2006 09:09:29 -0000	1.6
--- eigensystem.h	28 Oct 2006 00:30:44 -0000	1.7
***************
*** 7,10 ****
--- 7,11 ----
  #include <scitbx/array_family/accessors/c_grid.h>
  #include <scitbx/sym_mat3.h>
+ #include <scitbx/sym_mat2.h>
  #include <boost/scoped_array.hpp>
  
***************
*** 226,229 ****
--- 227,234 ----
          FloatType epsilon=1.e-10);
  
+       real_symmetric(
+         scitbx::sym_mat2<FloatType> const& m,
+         FloatType epsilon=1.e-10);
+ 
        //! The list of eigenvectors.
        af::versa<FloatType, af::c_grid<2> >
***************
*** 302,305 ****
--- 307,328 ----
    }
  
+   template <typename FloatType>
+   real_symmetric<FloatType>::real_symmetric(
+     scitbx::sym_mat2<FloatType> const& m,
+     FloatType epsilon)
+   {
+     FloatType a[3];
+     a[0] = m(0,0);
+     a[1] = m(1,0);
+     a[2] = m(1,1);
+     vectors_.resize(af::c_grid<2>(2, 2));
+     values_.resize(2);
+     detail::real_symmetric_given_lower_triangle(
+       a,
+       std::size_t(2),
+       vectors_.begin(),
+       values_.begin(),
+       epsilon);
+   }
  }}} // namespace scitbx::math::eigensystem
  

Index: principal_axes_of_inertia.h
===================================================================
RCS file: /cvsroot/cctbx/scitbx/include/scitbx/math/principal_axes_of_inertia.h,v
retrieving revision 1.8
retrieving revision 1.9
diff -C2 -d -r1.8 -r1.9
*** principal_axes_of_inertia.h	24 Oct 2006 18:32:04 -0000	1.8
--- principal_axes_of_inertia.h	28 Oct 2006 00:30:44 -0000	1.9
***************
*** 141,144 ****
--- 141,255 ----
    };
  
+   /*! Principal axes of inertia given discrete points in 2-dimensional
+       space and optionally weights.
+    */
+   template <typename FloatType=double>
+   class principal_axes_of_inertia_2d
+   {
+     private:
+       vec2<FloatType> center_of_mass_;
+       sym_mat2<FloatType> inertia_tensor_;
+       eigensystem::real_symmetric<FloatType> eigensystem_;
+ 
+     public:
+       //! Default constructor. Some data members are not initialized!
+       principal_axes_of_inertia_2d() {}
+ 
+       //! Intitialization given discrete points with unit weights.
+       principal_axes_of_inertia_2d(
+         af::const_ref<vec2<FloatType> > const& points)
+       :
+         center_of_mass_(0,0),
+         inertia_tensor_(0,0,0)
+       {
+         if (points.size() != 0) {
+           for(std::size_t i_p=0;i_p<points.size();i_p++) {
+             center_of_mass_ += points[i_p];
+           }
+           center_of_mass_ /= static_cast<FloatType>(points.size());
+           for(std::size_t i_p=0;i_p<points.size();i_p++) {
+             vec2<FloatType> p = points[i_p] - center_of_mass_;
+             vec2<FloatType> pp(p[0]*p[0], p[1]*p[1]);
+             inertia_tensor_(0,0) += pp[1];
+             inertia_tensor_(1,1) += pp[0];
+             inertia_tensor_(0,1) -= p[0] * p[1];
+           }
+         }
+         eigensystem_ = math::eigensystem::real_symmetric<FloatType>(
+           inertia_tensor_);
+       }
+ 
+       //! Intitialization given discrete points and weights.
+       principal_axes_of_inertia_2d(
+         af::const_ref<vec2<FloatType> > const& points,
+         af::const_ref<FloatType> const& weights)
+       :
+         center_of_mass_(0,0),
+         inertia_tensor_(0,0,0)
+       {
+         SCITBX_ASSERT(weights.size() == points.size());
+         FloatType sum_weights = 0;
+         for(std::size_t i_p=0;i_p<points.size();i_p++) {
+           FloatType w = weights[i_p];
+           if (w < 0) {
+             throw std::runtime_error(detail::report_negative_weight(
+               static_cast<double>(w), __FILE__, __LINE__));
+           }
+           center_of_mass_ += w * points[i_p];
+           sum_weights += w;
+         }
+         if (sum_weights != 0) {
+           center_of_mass_ /= sum_weights;
+           for(std::size_t i_p=0;i_p<points.size();i_p++) {
+             vec2<FloatType> p = points[i_p] - center_of_mass_;
+             vec2<FloatType> pp(p[0]*p[0], p[1]*p[1]);
+             FloatType w = weights[i_p];
+             inertia_tensor_(0,0) += w * pp[1];
+             inertia_tensor_(1,1) += w * pp[0];
+             inertia_tensor_(0,1) -= w * p[0] * p[1];
+ 
+           }
+         }
+         eigensystem_ = math::eigensystem::real_symmetric<FloatType>(
+           inertia_tensor_);
+       }
+ 
+       //! Center of mass.
+       /*! The weighted average of the coordinates of the points as passed
+           to the constructor.
+        */
+       vec2<FloatType> const&
+       center_of_mass() const { return center_of_mass_; }
+ 
+       //! Real-symmetric 2x2 inertia tensor.
+       /*! See e.g. http://kwon3d.com/theory/moi/iten.html or
+           search for "inertia tensor" at http://www.google.com/
+           or another search engine.
+        */
+       sym_mat2<FloatType> const&
+       inertia_tensor() const { return inertia_tensor_; }
+ 
+       //! Eigenvectors and eigenvalues of inertia_tensor().
+       math::eigensystem::real_symmetric<FloatType> const&
+       eigensystem() const { return eigensystem_; }
+ 
+       //! Distance from center of inertia tensor ellipsoid to its surface.
+       /*! unit_direction must be a vector of length 1. This is not checked
+           to minimize the runtime overhead.
+           If the inertia tensor is degenerate the result is 0.
+        */
+       FloatType
+       distance_to_inertia_ellipsoid_surface(
+         vec2<FloatType> const& unit_direction) const
+       {
+         FloatType d = inertia_tensor_.determinant();
+         if (d == static_cast<FloatType>(0)) return 0;
+         FloatType denom = (  inertia_tensor_.co_factor_matrix_transposed()
+                            * unit_direction).length();
+         if (denom == static_cast<FloatType>(0)) return 0;
+         return d / denom;
+       }
+   };
+ 
  }} // namespace scitbx::math