MolDrawing.h
Go to the documentation of this file.
// $Id$
//
// Copyright (C) 2009-2012 Greg Landrum
//
// @@ All Rights Reserved @@
// This file is part of the RDKit.
// The contents are covered by the terms of the BSD license
// which is included in the file license.txt, found at the root
// of the RDKit source tree.
//
// Includes contributions from Dave Cosgrove (davidacosgroveaz@gmail.com)
//
#ifndef _RD_MOLDRAWING_H_
#define _RD_MOLDRAWING_H_
#include <vector>
#include <boost/foreach.hpp>
#include <boost/lexical_cast.hpp>
#include <GraphMol/RDKitBase.h>
#include <GraphMol/Depictor/RDDepictor.h>
#include <Geometry/point.h>
/***********
Return Format: vector of ints
RESOLUTION dots_per_angstrom
BOUNDS x1 y1 x2 y2
LINE width dashed atom1_atnum atom2_atnum x1 y1 x2 y2
WEDGE dashed atom1_atnum atom2_atnum x1 y1 x2 y2 x3 y3
ATOM idx atnum x y num_chars char1-charx orient
*************/
namespace RDKit {
namespace Drawing {
typedef int ElementType;
typedef enum { LINE = 1, WEDGE, ATOM, BOUNDS, RESOLUTION } PrimType;
typedef enum { C = 0, N, E, S, W } OrientType;
namespace detail {
// **************************************************************************
void drawLine(std::vector<ElementType> &res, int atnum1, int atnum2,
int lineWidth, int dashed, double x1, double y1, double x2,
double y2) {
res.push_back(LINE);
res.push_back(static_cast<ElementType>(lineWidth));
res.push_back(dashed);
res.push_back(static_cast<ElementType>(atnum1));
res.push_back(static_cast<ElementType>(atnum2));
res.push_back(static_cast<ElementType>(x1));
res.push_back(static_cast<ElementType>(y1));
res.push_back(static_cast<ElementType>(x2));
res.push_back(static_cast<ElementType>(y2));
}
std::pair<std::string, OrientType> getAtomSymbolAndOrientation(
const Atom &atom, RDGeom::Point2D nbrSum) {
std::string symbol = "";
OrientType orient = C;
int isotope = atom.getIsotope();
if (atom.getAtomicNum() != 6 || atom.getFormalCharge() != 0 || isotope ||
atom.getNumRadicalElectrons() != 0 ||
atom.hasProp(common_properties::molAtomMapNumber) ||
atom.getDegree() == 0) {
symbol = atom.getSymbol();
bool leftToRight = true;
if (atom.getDegree() == 1 && nbrSum.x > 0) {
leftToRight = false;
}
if (isotope) {
symbol = boost::lexical_cast<std::string>(isotope) + symbol;
}
if (atom.hasProp(common_properties::molAtomMapNumber)) {
std::string mapNum;
atom.getProp(common_properties::molAtomMapNumber, mapNum);
symbol += ":" + mapNum;
}
int nHs = atom.getTotalNumHs();
if (nHs > 0) {
std::string h = "H";
if (nHs > 1) {
h += boost::lexical_cast<std::string>(nHs);
}
if (leftToRight)
symbol += h;
else
symbol = h + symbol;
}
if (atom.getFormalCharge() != 0) {
int chg = atom.getFormalCharge();
std::string sgn = "+";
if (chg < 0) {
sgn = "-";
}
chg = abs(chg);
if (chg > 1) {
sgn += boost::lexical_cast<std::string>(chg);
}
if (leftToRight)
symbol += sgn;
else
symbol = sgn + symbol;
}
if (atom.getDegree() == 1) {
double islope = 0;
if (fabs(nbrSum.y) > 1) {
islope = nbrSum.x / fabs(nbrSum.y);
} else {
islope = nbrSum.x;
}
if (fabs(islope) > .85) {
if (islope > 0) {
orient = W;
} else {
orient = E;
}
} else {
if (nbrSum.y > 0) {
orient = N;
} else {
orient = S;
}
}
}
}
return std::make_pair(symbol, orient);
}
} // end of detail namespace
// **************************************************************************
std::vector<ElementType> DrawMol(const ROMol &mol, int confId = -1,
const std::vector<int> *highlightAtoms = 0,
bool includeAtomCircles = false,
unsigned int dotsPerAngstrom = 100,
double dblBondOffset = 0.3,
double dblBondLengthFrac = 0.8,
double angstromsPerChar = 0.20) {
if (!mol.getRingInfo()->isInitialized()) {
MolOps::findSSSR(mol);
}
std::vector<ElementType> res;
res.push_back(RESOLUTION);
res.push_back(static_cast<ElementType>(dotsPerAngstrom));
const Conformer &conf = mol.getConformer(confId);
const RDGeom::POINT3D_VECT &locs = conf.getPositions();
// get atom symbols and orientations
// (we need them for the bounding box calculation)
std::vector<std::pair<std::string, OrientType> > atomSymbols;
ROMol::VERTEX_ITER bAts, eAts;
boost::tie(bAts, eAts) = mol.getVertices();
while (bAts != eAts) {
ROMol::OEDGE_ITER nbr, endNbrs;
RDGeom::Point2D nbrSum(0, 0);
boost::tie(nbr, endNbrs) = mol.getAtomBonds(mol[*bAts].get());
RDGeom::Point2D a1(locs[mol[*bAts]->getIdx()].x,
locs[mol[*bAts]->getIdx()].y);
while (nbr != endNbrs) {
const BOND_SPTR bond = mol[*nbr];
++nbr;
int a2Idx = bond->getOtherAtomIdx(mol[*bAts]->getIdx());
RDGeom::Point2D a2(locs[a2Idx].x, locs[a2Idx].y);
nbrSum += a2 - a1;
}
atomSymbols.push_back(
detail::getAtomSymbolAndOrientation(*mol[*bAts], nbrSum));
++bAts;
}
//------------
// do the bounding box
//------------
double minx = 1e6, miny = 1e6, maxx = -1e6, maxy = -1e6;
for (unsigned int i = 0; i < mol.getNumAtoms(); ++i) {
RDGeom::Point3D pt = locs[i];
std::string symbol;
OrientType orient;
boost::tie(symbol, orient) = atomSymbols[i];
if (symbol != "") {
// account for a possible expansion of the bounding box by the symbol
if (pt.x <= minx) {
switch (orient) {
case C:
case N:
case S:
case E:
minx = pt.x - symbol.size() / 2 * angstromsPerChar;
break;
case W:
minx = pt.x - symbol.size() * angstromsPerChar;
break;
}
}
if (pt.x >= maxx) {
switch (orient) {
case C:
case N:
case S:
case W:
maxx = pt.x + symbol.size() / 2 * angstromsPerChar;
break;
case E:
maxx = pt.x + symbol.size() * angstromsPerChar;
break;
}
}
if (pt.y <= miny) {
miny = pt.y - 1.5 * angstromsPerChar;
}
if (pt.y >= maxy) {
maxy = pt.y + angstromsPerChar;
}
} else {
minx = std::min(pt.x, minx);
miny = std::min(pt.y, miny);
maxx = std::max(pt.x, maxx);
maxy = std::max(pt.y, maxy);
}
}
double dimx = (maxx - minx), dimy = (maxy - miny);
res.push_back(BOUNDS);
res.push_back(static_cast<ElementType>(dotsPerAngstrom * 0));
res.push_back(static_cast<ElementType>(dotsPerAngstrom * 0));
res.push_back(static_cast<ElementType>(dotsPerAngstrom * dimx));
res.push_back(static_cast<ElementType>(dotsPerAngstrom * dimy));
// loop over atoms:
boost::tie(bAts, eAts) = mol.getVertices();
while (bAts != eAts) {
int a1Idx = mol[*bAts]->getIdx();
RDGeom::Point2D a1(locs[a1Idx].x - minx, locs[a1Idx].y - miny);
ROMol::OEDGE_ITER nbr, endNbrs;
RDGeom::Point2D nbrSum(0, 0);
boost::tie(nbr, endNbrs) = mol.getAtomBonds(mol[*bAts].get());
while (nbr != endNbrs) {
const BOND_SPTR bond = mol[*nbr];
++nbr;
int a2Idx = bond->getOtherAtomIdx(a1Idx);
int lineWidth = 1;
if (highlightAtoms &&
std::find(highlightAtoms->begin(), highlightAtoms->end(), a1Idx) !=
highlightAtoms->end() &&
std::find(highlightAtoms->begin(), highlightAtoms->end(), a2Idx) !=
highlightAtoms->end()) {
lineWidth = 3;
}
RDGeom::Point2D a2(locs[a2Idx].x - minx, locs[a2Idx].y - miny);
nbrSum += a2 - a1;
if (a2Idx < a1Idx) continue;
// draw bond from a1 to a2.
int atnum1 = mol[*bAts]->getAtomicNum();
int atnum2 = mol.getAtomWithIdx(a2Idx)->getAtomicNum();
if (!mol.getRingInfo()->numBondRings(bond->getIdx()) &&
bond->getBondType() != Bond::AROMATIC) {
// acyclic bonds
RDGeom::Point2D obv = a2 - a1;
RDGeom::Point2D perp = obv;
perp.rotate90();
perp.normalize();
if (bond->getBondType() == Bond::DOUBLE ||
bond->getBondType() == Bond::TRIPLE) {
RDGeom::Point2D startP = a1, endP = a2;
if (bond->getBondType() == Bond::TRIPLE) {
perp *= dblBondOffset;
startP += (obv * (1. - dblBondLengthFrac) / 2);
endP -= (obv * (1. - dblBondLengthFrac) / 2);
} else {
perp *= 0.5 * dblBondOffset;
}
detail::drawLine(res, atnum1, atnum2, lineWidth, 0,
dotsPerAngstrom * (startP.x + perp.x),
dotsPerAngstrom * (startP.y + perp.y),
dotsPerAngstrom * (endP.x + perp.x),
dotsPerAngstrom * (endP.y + perp.y));
if (bond->getBondType() != Bond::AROMATIC) {
detail::drawLine(res, atnum1, atnum2, lineWidth, 0,
dotsPerAngstrom * (startP.x - perp.x),
dotsPerAngstrom * (startP.y - perp.y),
dotsPerAngstrom * (endP.x - perp.x),
dotsPerAngstrom * (endP.y - perp.y));
} else {
detail::drawLine(res, atnum1, atnum2, lineWidth, 1,
dotsPerAngstrom * (startP.x - perp.x),
dotsPerAngstrom * (startP.y - perp.y),
dotsPerAngstrom * (endP.x - perp.x),
dotsPerAngstrom * (endP.y - perp.y));
}
}
if (bond->getBondType() == Bond::SINGLE ||
bond->getBondType() == Bond::TRIPLE) {
detail::drawLine(res, atnum1, atnum2, lineWidth, 0,
dotsPerAngstrom * (a1.x), dotsPerAngstrom * (a1.y),
dotsPerAngstrom * (a2.x), dotsPerAngstrom * (a2.y));
} else if (bond->getBondType() != Bond::DOUBLE) {
detail::drawLine(res, atnum1, atnum2, lineWidth, 2,
dotsPerAngstrom * (a1.x), dotsPerAngstrom * (a1.y),
dotsPerAngstrom * (a2.x), dotsPerAngstrom * (a2.y));
}
} else {
// cyclic bonds
detail::drawLine(res, atnum1, atnum2, lineWidth, 0,
dotsPerAngstrom * a1.x, dotsPerAngstrom * a1.y,
dotsPerAngstrom * a2.x, dotsPerAngstrom * a2.y);
if (bond->getBondType() == Bond::DOUBLE ||
bond->getBondType() == Bond::AROMATIC ||
bond->getBondType() == Bond::TRIPLE) {
RDGeom::Point2D obv = a2 - a1;
RDGeom::Point2D perp = obv;
perp.rotate90();
perp.normalize();
if ((bond->getBondType() == Bond::DOUBLE ||
bond->getBondType() == Bond::AROMATIC) &&
mol.getRingInfo()->numBondRings(bond->getIdx())) {
// we're in a ring... we might need to flip sides:
ROMol::OEDGE_ITER nbr2, endNbrs2;
boost::tie(nbr2, endNbrs2) = mol.getAtomBonds(mol[*bAts].get());
while (nbr2 != endNbrs2) {
const BOND_SPTR bond2 = mol[*nbr2];
++nbr2;
if (bond2->getIdx() == bond->getIdx() ||
!mol.getRingInfo()->numBondRings(bond2->getIdx()))
continue;
bool sharedRing = false;
BOOST_FOREACH (const INT_VECT &ring,
mol.getRingInfo()->bondRings()) {
if (std::find(ring.begin(), ring.end(), bond->getIdx()) !=
ring.end() &&
std::find(ring.begin(), ring.end(), bond2->getIdx()) !=
ring.end()) {
sharedRing = true;
break;
}
}
if (sharedRing) {
// these two bonds share a ring.
int a3Idx = bond2->getOtherAtomIdx(a1Idx);
if (a3Idx != a2Idx) {
RDGeom::Point2D a3(locs[a3Idx].x - minx,
locs[a3Idx].y - miny);
RDGeom::Point2D obv2 = a3 - a1;
if (obv2.dotProduct(perp) < 0) {
perp *= -1;
}
}
}
}
}
perp *= dblBondOffset;
RDGeom::Point2D offsetStart =
a1 + obv * (.5 * (1. - dblBondLengthFrac));
obv *= dblBondLengthFrac;
detail::drawLine(res, atnum1, atnum2, lineWidth,
(bond->getBondType() == Bond::AROMATIC),
dotsPerAngstrom * (offsetStart.x + perp.x),
dotsPerAngstrom * (offsetStart.y + perp.y),
dotsPerAngstrom * (offsetStart.x + obv.x + perp.x),
dotsPerAngstrom * (offsetStart.y + obv.y + perp.y));
}
}
}
std::string symbol;
OrientType orient;
boost::tie(symbol, orient) = atomSymbols[a1Idx];
if (symbol != "" || includeAtomCircles) {
res.push_back(ATOM);
res.push_back(mol[*bAts]->getAtomicNum());
res.push_back(static_cast<ElementType>(dotsPerAngstrom * a1.x));
res.push_back(static_cast<ElementType>(dotsPerAngstrom * a1.y));
res.push_back(static_cast<ElementType>(symbol.length()));
if (symbol.length()) {
BOOST_FOREACH (char c, symbol) {
res.push_back(static_cast<ElementType>(c));
}
}
res.push_back(static_cast<ElementType>(orient));
}
++bAts;
}
return res;
}
std::vector<int> MolToDrawing(const RDKit::ROMol &mol,
const std::vector<int> *highlightAtoms = 0,
bool kekulize = true,
bool includeAtomCircles = false) {
RDKit::RWMol *cp = new RDKit::RWMol(mol);
if (kekulize) {
try {
RDKit::MolOps::Kekulize(*cp);
} catch (...) {
delete cp;
cp = new RDKit::RWMol(mol);
}
}
if (!mol.getNumConformers()) {
RDDepict::compute2DCoords(*cp);
}
std::vector<int> drawing =
DrawMol(*cp, -1, highlightAtoms, includeAtomCircles);
delete cp;
return drawing;
}
} // end of namespace Drawing
} // end of namespace RDKit
#endif
RDKit::BOND_SPTR
boost::shared_ptr< Bond > BOND_SPTR
Definition: ROMol.h:41
RDKit::ROMol::getConformer
const Conformer & getConformer(int id=-1) const
RDKit::Drawing::N
Definition: MolDrawing.h:41
RDKit::Drawing::detail::getAtomSymbolAndOrientation
std::pair< std::string, OrientType > getAtomSymbolAndOrientation(const Atom &atom, RDGeom::Point2D nbrSum)
Definition: MolDrawing.h:58
RDKit::MolOps::Kekulize
void Kekulize(RWMol &mol, bool markAtomsBonds=true, unsigned int maxBackTracks=100)
Kekulizes the molecule.
RDKit::Drawing::MolToDrawing
std::vector< int > MolToDrawing(const RDKit::ROMol &mol, const std::vector< int > *highlightAtoms=0, bool kekulize=true, bool includeAtomCircles=false)
Definition: MolDrawing.h:396
RDKit::ROMol::getVertices
ATOM_ITER_PAIR getVertices()
returns an iterator pair for looping over all Atoms
RDKit::Drawing::WEDGE
Definition: MolDrawing.h:40
RDKit::RingInfo::numBondRings
unsigned int numBondRings(unsigned int idx) const
returns the number of rings bond idx is involved in
RDKit::MolOps::findSSSR
int findSSSR(const ROMol &mol, std::vector< std::vector< int > > &res)
finds a molecule&#39;s Smallest Set of Smallest Rings
RDGeom::Point2D::rotate90
void rotate90()
Definition: point.h:330
RDGeom::Point2D::y
double y
Definition: point.h:256
RDKit::RWMol
RWMol is a molecule class that is intended to be edited.
Definition: RWMol.h:30
RDKit::ROMol::getNumAtoms
unsigned int getNumAtoms(bool onlyExplicit=1) const
returns our number of atoms
RDKit::RingInfo::bondRings
const VECT_INT_VECT & bondRings() const
returns our bond-rings vectors
Definition: RingInfo.h:126
point.h
RDKit::common_properties::molAtomMapNumber
const std::string molAtomMapNumber
RDKit::Atom::getTotalNumHs
unsigned int getTotalNumHs(bool includeNeighbors=false) const
returns the total number of Hs (implicit and explicit) that this Atom is bound to ...
RDKit::Drawing::C
Definition: MolDrawing.h:41
RDGeom::Point2D
Definition: point.h:254
RDKit::ROMol::getNumConformers
unsigned int getNumConformers() const
Definition: ROMol.h:363
RDKit::Atom::getNumRadicalElectrons
unsigned int getNumRadicalElectrons() const
returns the number of radical electrons for this Atom
Definition: Atom.h:193
RDKit::Drawing::S
Definition: MolDrawing.h:41
RDGeom::POINT3D_VECT
std::vector< Point3D > POINT3D_VECT
Definition: point.h:507
RDKitBase.h
pulls in the core RDKit functionality
RDKit::ROMol
ROMol is a molecule class that is intended to have a fixed topology.
Definition: ROMol.h:102
RDKit::Drawing::RESOLUTION
Definition: MolDrawing.h:40
RDKit::Atom::hasProp
bool hasProp(const char *key) const
returns whether or not we have a property with name key
Definition: Atom.h:396
RDKit::ROMol::getRingInfo
RingInfo * getRingInfo() const
Definition: ROMol.h:374
RDKit::Conformer::getPositions
const RDGeom::POINT3D_VECT & getPositions() const
Get a const reference to the vector of atom positions.
RDKit::Atom::getIsotope
unsigned int getIsotope() const
returns our isotope number
Definition: Atom.h:225
RDKit::Bond::DOUBLE
Definition: Bond.h:59
RDGeom::Point2D::dotProduct
double dotProduct(const Point2D &other) const
Definition: point.h:347
detail
Definition: types.h:23
RDKit::Drawing::LINE
Definition: MolDrawing.h:40
RDDepictor.h
RDKit::INT_VECT
std::vector< int > INT_VECT
Definition: types.h:146
RDGeom::Point2D::normalize
void normalize()
Definition: point.h:324
RDKit::Drawing::BOUNDS
Definition: MolDrawing.h:40
RDKit::Drawing::E
Definition: MolDrawing.h:41
RDKit::Atom::getFormalCharge
int getFormalCharge() const
returns the formal charge of this atom
Definition: Atom.h:199
RDKit
Includes a bunch of functionality for handling Atom and Bond queries.
Definition: Atom.h:28
RDKit::Drawing::ATOM
Definition: MolDrawing.h:40
RDKit::Drawing::DrawMol
std::vector< ElementType > DrawMol(const ROMol &mol, int confId=-1, const std::vector< int > *highlightAtoms=0, bool includeAtomCircles=false, unsigned int dotsPerAngstrom=100, double dblBondOffset=0.3, double dblBondLengthFrac=0.8, double angstromsPerChar=0.20)
Definition: MolDrawing.h:133
RDKit::Drawing::OrientType
OrientType
Definition: MolDrawing.h:41
RDKit::Atom::getAtomicNum
int getAtomicNum() const
returns our atomic number
Definition: Atom.h:115
RDKit::ROMol::getAtomBonds
OBOND_ITER_PAIR getAtomBonds(Atom const *at) const
provides access to all Bond objects connected to an Atom
RDGeom::Point3D::y
double y
Definition: point.h:47
RDGeom::Point3D
Definition: point.h:45
RDKit::Drawing::W
Definition: MolDrawing.h:41
RDKit::Bond::TRIPLE
Definition: Bond.h:60
RDDepict::compute2DCoords
unsigned int compute2DCoords(RDKit::ROMol &mol, const RDGeom::INT_POINT2D_MAP *coordMap=0, bool canonOrient=false, bool clearConfs=true, unsigned int nFlipsPerSample=0, unsigned int nSamples=0, int sampleSeed=0, bool permuteDeg4Nodes=false)
Generate 2D coordinates (a depiction) for a molecule.
RDKit::Conformer
The class for representing 2D or 3D conformation of a molecule.
Definition: Conformer.h:41
RDKit::Drawing::PrimType
PrimType
Definition: MolDrawing.h:40
RDKit::Bond::AROMATIC
Definition: Bond.h:69
RDGeom::Point2D::x
double x
Definition: point.h:256
RDKit::Atom::getDegree
unsigned int getDegree() const
RDKit::ROMol::getAtomWithIdx
Atom * getAtomWithIdx(unsigned int idx)
returns a pointer to a particular Atom
RDGeom::Point3D::x
double x
Definition: point.h:47
RDKit::Atom::getSymbol
std::string getSymbol() const
returns our symbol (determined by our atomic number)
RDKit::Drawing::ElementType
int ElementType
Definition: MolDrawing.h:38
RDKit::Bond::SINGLE
Definition: Bond.h:58
RDKit::Atom
The class for representing atoms.
Definition: Atom.h:67
RDKit::Atom::getProp
void getProp(const char *key, T &res) const
allows retrieval of a particular property value
Definition: Atom.h:363
RDKit::RingInfo::isInitialized
bool isInitialized() const
checks to see if we&#39;ve been properly initialized
Definition: RingInfo.h:39
RDKit::Drawing::detail::drawLine
void drawLine(std::vector< ElementType > &res, int atnum1, int atnum2, int lineWidth, int dashed, double x1, double y1, double x2, double y2)
Definition: MolDrawing.h:45