kalzium/libavogadro-kalzium/src/primitive.cpp

   1 /**********************************************************************
   2   Primitive - Wrapper class around the OpenBabel classes
   3
   4   Copyright (C) 2007 Donald Ephraim Curtis
   5
   6   This file is part of the Avogadro molecular editor project.
   7   For more information, see <http://avogadro.sourceforge.net/>
   8
   9   Avogadro is free software; you can redistribute it and/or modify
  10   it under the terms of the GNU General Public License as published by
  11   the Free Software Foundation; either version 2 of the License, or
  12   (at your option) any later version.
  13
  14   Avogadro is distributed in the hope that it will be useful,
  15   but WITHOUT ANY WARRANTY; without even the implied warranty of
  16   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17   GNU General Public License for more details.
  18
  19   You should have received a copy of the GNU General Public License
  20   along with this program; if not, write to the Free Software
  21   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  22   02110-1301, USA.
  23  **********************************************************************/
  24
  25 #include <config.h>
  26
  27 #include <avogadro/primitive.h>
  28 #include <eigen/regression.h>
  29
  30 #include <QReadWriteLock>
  31
  32 using namespace OpenBabel;
  33
  34 namespace Avogadro {
  35
  36   class PrimitivePrivate {
  37     public:
  38       PrimitivePrivate() : type(Primitive::OtherType) {};
  39
  40       enum Primitive::Type type;
  41       QReadWriteLock lock;
  42   };
  43
  44   Primitive::Primitive(QObject *parent) : QObject(parent), d_ptr(new PrimitivePrivate) {}
  45
  46   Primitive::Primitive(enum Type type, QObject *parent) : QObject(parent), d_ptr(new PrimitivePrivate)
  47   {
  48     Q_D(Primitive);
  49     d->type = type;
  50   }
  51
  52   Primitive::Primitive(PrimitivePrivate &dd, QObject *parent) : QObject(parent), d_ptr(&dd) {}
  53
  54   Primitive::Primitive(PrimitivePrivate &dd, enum Type type, QObject *parent) : QObject(parent), d_ptr(&dd)
  55   {
  56     Q_D(Primitive);
  57     d->type = type;
  58   }
  59
  60   Primitive::~Primitive()
  61   {
  62     delete d_ptr;
  63   }
  64
  65   enum Primitive::Type Primitive::type() const
  66   {
  67     Q_D(const Primitive);
  68     return d->type;
  69   }
  70
  71   QReadWriteLock *Primitive::lock()
  72   {
  73     Q_D(Primitive);
  74     return &d->lock;
  75   }
  76
  77   void Primitive::update()
  78   {
  79     emit updated();
  80   }
  81
  82   class MoleculePrivate : public PrimitivePrivate {
  83     public:
  84       MoleculePrivate() : PrimitivePrivate(), farthestAtom(0), invalidGeomInfo(true) {}
  85       mutable Eigen::Vector3d       center;
  86       mutable Eigen::Vector3d       normalVector;
  87       mutable double                radius;
  88       mutable Atom *                farthestAtom;
  89       mutable bool                  invalidGeomInfo;
  90   };
  91
  92   Molecule::Molecule(QObject *parent) : Primitive(*new MoleculePrivate, MoleculeType, parent), OpenBabel::OBMol()
  93   {
  94     connect(this, SIGNAL(updated()), this, SLOT(updatePrimitive()));
  95   }
  96
  97   Molecule::Molecule(const Molecule &other) : Primitive(*new MoleculePrivate, MoleculeType, other.parent()), OpenBabel::OBMol(other)
  98   {
  99     connect(this, SIGNAL(updated()), this, SLOT(updatePrimitive()));
 100   }
 101
 102   Molecule::~Molecule()
 103   {
 104   }
 105
 106   Atom * Molecule::CreateAtom()
 107   {
 108     Atom *atom = new Atom(this);
 109     connect(atom, SIGNAL(updated()), this, SLOT(updatePrimitive()));
 110     emit primitiveAdded(atom);
 111     return(atom);
 112   }
 113
 114   Bond * Molecule::CreateBond()
 115   {
 116     Bond *bond = new Bond(this);
 117     connect(bond, SIGNAL(updated()), this, SLOT(updatePrimitive()));
 118     emit primitiveAdded(bond);
 119     return(bond);
 120   }
 121
 122   Residue * Molecule::CreateResidue()
 123   {
 124     Residue *residue = new Residue(this);
 125     connect(residue, SIGNAL(updated()), this, SLOT(updatePrimitive()));
 126     emit primitiveAdded(residue);
 127     return(residue);
 128   }
 129
 130   void Molecule::DestroyAtom(OpenBabel::OBAtom *obatom)
 131   {
 132     Atom *atom = static_cast<Atom *>(obatom);
 133     if(atom) {
 134       emit primitiveRemoved(atom);
 135       atom->deleteLater();
 136     }
 137   }
 138
 139   void Molecule::DestroyBond(OpenBabel::OBBond *obbond)
 140   {
 141     Bond *bond = static_cast<Bond *>(obbond);
 142     if(bond) {
 143       emit primitiveRemoved(bond);
 144       bond->deleteLater();
 145     }
 146   }
 147
 148   void Molecule::DestroyResidue(OpenBabel::OBResidue *obresidue)
 149   {
 150     Residue *residue = static_cast<Residue *>(obresidue);
 151     if(residue) {
 152       emit primitiveRemoved(residue);
 153       residue->deleteLater();
 154     }
 155   }
 156
 157   void Molecule::updatePrimitive()
 158   {
 159     Q_D(Molecule);
 160     Primitive *primitive = qobject_cast<Primitive *>(sender());
 161     d->invalidGeomInfo = true;
 162     emit primitiveUpdated(primitive);
 163   }
 164
 165   void Molecule::update()
 166   {
 167     Q_D(Molecule);
 168     d->invalidGeomInfo = true;
 169     emit updated();
 170   }
 171
 172   const Eigen::Vector3d & Molecule::center() const
 173   {
 174     Q_D(const Molecule);
 175     if( d->invalidGeomInfo ) computeGeomInfo();
 176     return d->center;
 177   }
 178
 179   const Eigen::Vector3d & Molecule::normalVector() const
 180   {
 181     Q_D(const Molecule);
 182     if( d->invalidGeomInfo ) computeGeomInfo();
 183     return d->normalVector;
 184   }
 185
 186   const double & Molecule::radius() const
 187   {
 188     Q_D(const Molecule);
 189     if( d->invalidGeomInfo ) computeGeomInfo();
 190     return d->radius;
 191   }
 192
 193   const Atom * Molecule::farthestAtom() const
 194   {
 195     Q_D(const Molecule);
 196     if( d->invalidGeomInfo ) computeGeomInfo();
 197     return d->farthestAtom;
 198   }
 199
 200   Molecule &Molecule::operator=(const Molecule& other)
 201   {
 202     OpenBabel::OBMol::operator=(other);
 203         return *this;
 204   }
 205
 206   Molecule &Molecule::operator+=(const Molecule& other)
 207   {
 208     OpenBabel::OBMol::operator+=(other);
 209
 210     return *this;
 211   }
 212
 213   void Molecule::computeGeomInfo() const
 214   {
 215 //     MoleculePrivate *d = reinterpret_cast<MoleculePrivate *>(d_ptr);
 216     Q_D(const Molecule);
 217     d->invalidGeomInfo = true;
 218     d->farthestAtom = 0;
 219     d->center.loadZero();
 220     d->normalVector.loadZero();
 221     d->radius = 0.0;
 222     if( NumAtoms() != 0 )
 223     {
 224       // compute center
 225       std::vector< OpenBabel::OBAtom * >::iterator atom_iterator;
 226       for( Atom* atom = (Atom*) const_cast<Molecule*>(this)->BeginAtom(atom_iterator); atom; atom = (Atom *) const_cast<Molecule*>(this)->NextAtom(atom_iterator) )
 227       {
 228         d->center += atom->pos();
 229       }
 230       d->center /= NumAtoms();
 231
 232       // compute the normal vector to the molecule's best-fitting plane
 233       Eigen::Vector3d * atomPositions = new Eigen::Vector3d[NumAtoms()];
 234       int i = 0;
 235       for( Atom* atom = (Atom*) const_cast<Molecule*>(this)->BeginAtom(atom_iterator); atom; atom = (Atom *) const_cast<Molecule*>(this)->NextAtom(atom_iterator) )
 236       {
 237         atomPositions[i++] = atom->pos();
 238       }
 239       Eigen::Vector4d planeCoeffs;
 240       Eigen::computeFittingHyperplane( NumAtoms(), atomPositions, &planeCoeffs );
 241       delete[] atomPositions;
 242       d->normalVector = Eigen::Vector3d( planeCoeffs.x(), planeCoeffs.y(), planeCoeffs.z() );
 243       d->normalVector.normalize();
 244
 245       // compute radius and the farthest atom
 246       d->radius = -1.0; // so that ( squaredDistanceToCenter > d->radius ) is true for at least one atom.
 247       for( Atom* atom = (Atom*) const_cast<Molecule*>(this)->BeginAtom(atom_iterator); atom; atom = (Atom *) const_cast<Molecule*>(this)->NextAtom(atom_iterator) )
 248       {
 249         double distanceToCenter = (atom->pos() - d->center).norm();
 250         if( distanceToCenter > d->radius )
 251         {
 252           d->radius = distanceToCenter;
 253           d->farthestAtom = atom;
 254         }
 255       }
 256     }
 257     d->invalidGeomInfo = false;
 258   }
 259
 260 }
 261
 262 #include "primitive.moc"