bounding_box.cc

//------------------------------------------------------------------------------
// This file is part of the OpenStructure project <www.openstructure.org>
//
// Copyright (C) 2008-2015 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
//------------------------------------------------------------------------------

/*
  Author: Marco Biasini
 */

#include <limits>
#include <Eigen/QR>
#include <ost/mol/mol.hh>
#include <ost/mol/bounding_box.hh>
namespace ost { namespace mol {

namespace {

typedef Eigen::Matrix3d EMat3;

class PrincipalAxisCalc : public EntityVisitor {
public:
  PrincipalAxisCalc(const EntityHandle& ent): 
    natoms_(ent.GetAtomCount())
  {
    mean_=ent.GetBounds().GetCenter();
    cov_.setZero();
  }
  PrincipalAxisCalc(const AtomHandleList& atoms):
    natoms_(atoms.size())
  { 
    cov_.setZero();    
    mean_=geom::Vec3();
    if (natoms_>0) {
      for (AtomHandleList::const_iterator i=atoms.begin(),
           e=atoms.end(); i!=e; ++i) {
        mean_+=i->GetPos();       
      }      
    }
    mean_/=natoms_;
  }
  PrincipalAxisCalc(const EntityView& ent): 
    natoms_(ent.GetAtomCount())
  {
    mean_=ent.GetBounds().GetCenter();
    cov_.setZero();
  }  
  virtual bool VisitAtom(const AtomHandle& atom)
  {
    this->UpdateCov(atom);
    return false;
  } 
  void UpdateCov(const AtomHandle& atom)
  {
    geom::Vec3 p=atom.GetPos();
    for (size_t i=0; i<3; ++i) {
      for (size_t j=0; j<3; ++j) {