pub_soft/rdkit
3
0
Fork 0
mirror of https://github.com/rdkit/rdkit.git synced 2026-06-03 21:44:30 +08:00
Code Issues Packages Projects Releases Wiki Activity
Files
master
rdkit/Code/GraphMol/MolDraw2D/MolDraw2DUtils.cpp
Paolo Tosco 805294c27f - move MolDrawOptions JSON parsing code to a separate compilation unit (#8783) 
- implement drawingExtentsInclude MolDrawOption taking as parameter DrawElement enum values
- add unit tests for drawingExtentsInclude option (C++, Python, CFFI)

Co-authored-by: ptosco <paolo.tosco@novartis.com>
2025-09-23 09:27:13 +02:00

417 lines
15 KiB
C++
Raw Permalink Blame History

//
// Copyright (C) 2016-2021 Greg Landrum and other RDKit contributors
//
// @@ 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.
//
#include <GraphMol/MolDraw2D/MolDraw2D.h>
#include <GraphMol/FileParsers/MolFileStereochem.h>
#include <GraphMol/MolTransforms/MolTransforms.h>
#include <GraphMol/FileParsers/FileParsers.h>
#include <GraphMol/Chirality.h>
#ifdef RDK_BUILD_CAIRO_SUPPORT
#include <GraphMol/MolDraw2D/MolDraw2DCairo.h>
#endif
#include <GraphMol/MolDraw2D/MolDraw2DSVG.h>
#include <GraphMol/MolDraw2D/MolDraw2DUtils.h>
#include <GraphMol/RWMol.h>
#include <GraphMol/MolOps.h>
#include <GraphMol/Depictor/RDDepictor.h>
#include <RDGeneral/BoostStartInclude.h>
#include <boost/lexical_cast.hpp>
#include <boost/property_tree/ptree.hpp>
#include <boost/property_tree/json_parser.hpp>
#include <boost/algorithm/string.hpp>
#include <RDGeneral/BoostEndInclude.h>
#include <RDGeneral/JSONHelpers.h>
#include <limits>
#include <cmath>
#include <sys/stat.h>
#include <Numerics/Conrec.h>
namespace RDKit {
namespace MolDraw2DUtils {
namespace {
bool isAtomCandForChiralH(const RWMol &mol, const Atom *atom) {
// conditions for needing a chiral H:
// - stereochem specified
// - in at least two rings
return mol.getRingInfo()->isInitialized() &&
mol.getRingInfo()->numAtomRings(atom->getIdx()) > 1u &&
(atom->getChiralTag() == Atom::CHI_TETRAHEDRAL_CCW ||
atom->getChiralTag() == Atom::CHI_TETRAHEDRAL_CW);
}
} // end of anonymous namespace
void prepareMolForDrawing(RWMol &mol, bool kekulize, bool addChiralHs,
bool wedgeBonds, bool forceCoords, bool wavyBonds) {
if (kekulize) {
RDLog::LogStateSetter blocker;
MolOps::KekulizeIfPossible(
mol, false); // kekulize, but keep the aromatic flags!
}
if (addChiralHs) {
std::vector<unsigned int> chiralAts;
for (auto atom : mol.atoms()) {
if (isAtomCandForChiralH(mol, atom)) {
chiralAts.push_back(atom->getIdx());
}
}
if (chiralAts.size()) {
bool addCoords = false;
if (!forceCoords && mol.getNumConformers()) {
addCoords = true;
}
MolOps::addHs(mol, false, addCoords, &chiralAts);
}
}
if (forceCoords || !mol.getNumConformers()) {
// compute 2D coordinates in a standard orientation:
const bool canonOrient = true;
RDDepict::compute2DCoords(mol, nullptr, canonOrient);
}
if (wedgeBonds) {
Chirality::wedgeMolBonds(mol, &mol.getConformer());
}
if (wavyBonds) {
addWavyBondsForStereoAny(mol);
}
}
void prepareAndDrawMolecule(MolDraw2D &drawer, const ROMol &mol,
const std::string &legend,
const std::vector<int> *highlight_atoms,
const std::vector<int> *highlight_bonds,
const std::map<int, DrawColour> *highlight_atom_map,
const std::map<int, DrawColour> *highlight_bond_map,
const std::map<int, double> *highlight_radii,
int confId, bool kekulize, bool addChiralHs,
bool wedgeBonds, bool forceCoords, bool wavyBonds) {
RWMol cpy(mol);
prepareMolForDrawing(cpy, kekulize, addChiralHs, wedgeBonds, forceCoords,
wavyBonds);
// having done the prepare, we don't want to do it again in drawMolecule.
bool old_prep_mol = drawer.drawOptions().prepareMolsBeforeDrawing;
drawer.drawOptions().prepareMolsBeforeDrawing = false;
drawer.drawMolecule(cpy, legend, highlight_atoms, highlight_bonds,
highlight_atom_map, highlight_bond_map, highlight_radii,
confId);
drawer.drawOptions().prepareMolsBeforeDrawing = old_prep_mol;
}
void contourAndDrawGrid(MolDraw2D &drawer, const double *grid,
const std::vector<double> &xcoords,
const std::vector<double> &ycoords, size_t nContours,
std::vector<double> &levels,
const ContourParams &params, const ROMol *mol) {
PRECONDITION(grid, "no data");
PRECONDITION(params.colourMap.size() > 1,
"colourMap must have at least two entries");
if (params.setScale) {
Point2D minP = {xcoords[0], ycoords[0]};
Point2D maxP = {xcoords.back(), ycoords.back()};
drawer.setScale(drawer.width(), drawer.height(), minP, maxP, mol);
}
size_t nX = xcoords.size();
size_t nY = ycoords.size();
double minV = std::numeric_limits<double>::max();
double maxV = std::numeric_limits<double>::lowest();
if (!levels.size() || params.fillGrid) {
for (size_t i = 0; i < nX; ++i) {
for (size_t j = 0; j < nY; ++j) {
minV = std::min(minV, grid[i * nY + j]);
maxV = std::max(maxV, grid[i * nY + j]);
}
}
if (!levels.size()) {
levels.resize(nContours);
for (size_t i = 0; i < nContours; ++i) {
levels[i] = minV + i * (maxV - minV) / (nContours - 1);
}
}
}
if (maxV <= minV) {
return;
}
const auto olw = drawer.lineWidth();
const auto odash = drawer.dash();
const auto ocolor = drawer.colour();
const auto ofill = drawer.fillPolys();
const auto owidth = drawer.lineWidth();
if (params.fillGrid) {
drawer.setFillPolys(true);
drawer.setLineWidth(1);
auto delta = (maxV - minV);
if (params.colourMap.size() > 2) {
// need to find how fractionally far we are from zero, not the min
if (-minV > maxV) {
delta = -minV;
} else {
delta = maxV;
}
}
for (size_t i = 0; i < nX - 1; ++i) {
for (size_t j = 0; j < nY - 1; ++j) {
auto gridV = grid[i * nY + j];
auto threshTest = gridV;
if (params.fillThresholdIsFraction) {
threshTest /= delta;
}
if (params.useFillThreshold &&
fabs(threshTest) < params.fillThreshold) {
continue;
}
auto fracV = (gridV - minV) / delta;
if (params.colourMap.size() > 2) {
// need to find how fractionally far we are from zero, not the min
fracV = gridV / delta;
if (fracV < 0) {
fracV *= -1;
}
}
auto c1 = (gridV < 0 || params.colourMap.size() == 2)
? params.colourMap[1]
: params.colourMap[1];
auto c2 = (gridV < 0 || params.colourMap.size() == 2)
? params.colourMap[0]
: params.colourMap[2];
auto c = c1 + (c2 - c1) * fracV;
// don't bother drawing boxes that are the same as the background color:
double tol = 0.01;
if (c.feq(drawer.drawOptions().backgroundColour, tol)) {
continue;
}
drawer.setColour(c);
Point2D p1 = {xcoords[i], ycoords[j]};
Point2D p2 = {xcoords[i + 1], ycoords[j + 1]};
drawer.drawRect(p1, p2);
}
}
}
if (nContours) {
if (nContours > levels.size()) {
throw ValueErrorException(
"nContours larger than the size of the level list");
}
std::vector<conrec::ConrecSegment> segs;
conrec::Contour(grid, 0, nX - 1, 0, nY - 1, xcoords.data(), ycoords.data(),
nContours, levels.data(), segs);
static DashPattern negDash{2., 2.};
static DashPattern posDash;
drawer.setColour(params.contourColour);
drawer.setLineWidth(params.contourWidth);
if (!params.drawAsLines) {
for (const auto &seg : segs) {
if (params.dashNegative && seg.isoVal < 0) {
drawer.setDash(negDash);
} else {
drawer.setDash(posDash);
}
drawer.drawLine(seg.p1, seg.p2);
}
} else {
drawer.setFillPolys(false);
auto lines =
conrec::connectLineSegments(segs, params.coordScaleForQuantization,
params.isovalScaleForQuantization);
for (const auto &pr : lines) {
auto [contour, val] = pr;
if (params.dashNegative && val < 0) {
drawer.setDash(negDash);
} else {
drawer.setDash(posDash);
}
if (contour.size() > 2) {
drawer.drawPolygon(contour);
} else if (contour.size() == 2) {
drawer.drawLine(contour[0], contour[1]);
}
}
}
drawer.setDash(odash);
drawer.setLineWidth(olw);
drawer.setColour(ocolor);
drawer.setFillPolys(ofill);
drawer.setLineWidth(owidth);
}
};
void contourAndDrawGaussians(MolDraw2D &drawer,
const std::vector<Point2D> &locs,
const std::vector<double> &weights,
const std::vector<double> &widths,
size_t nContours, std::vector<double> &levels,
const ContourParams &params, const ROMol *mol) {
PRECONDITION(locs.size() == weights.size(), "size mismatch");
PRECONDITION(locs.size() == widths.size(), "size mismatch");
// start by setting up the grid
if (params.setScale) {
Point2D minP, maxP;
minP.x = minP.y = std::numeric_limits<double>::max();
maxP.x = maxP.y = std::numeric_limits<double>::lowest();
for (const auto &loc : locs) {
minP.x = std::min(loc.x, minP.x);
minP.y = std::min(loc.y, minP.y);
maxP.x = std::max(loc.x, maxP.x);
maxP.y = std::max(loc.y, maxP.y);
}
Point2D dims = maxP - minP;
// Here, the drawOptions().padding is just used to extend the grid
// beyond the molecule. The actual padding round the image is added
// later.
minP.x -= drawer.drawOptions().padding * dims.x;
minP.y -= drawer.drawOptions().padding * dims.y;
maxP.x += drawer.drawOptions().padding * dims.x;
maxP.y += drawer.drawOptions().padding * dims.y;
if (params.extraGridPadding > 0) {
minP.x -= params.extraGridPadding;
minP.y -= params.extraGridPadding;
maxP.x += params.extraGridPadding;
maxP.y += params.extraGridPadding;
}
drawer.setScale(drawer.width(), drawer.height(), minP, maxP, mol);
}
size_t nx = (size_t)ceil(drawer.range().x / params.gridResolution) + 1;
size_t ny = (size_t)ceil(drawer.range().y / params.gridResolution) + 1;
std::vector<double> xcoords(nx);
for (size_t i = 0; i < nx; ++i) {
xcoords[i] = drawer.minPt().x + i * params.gridResolution;
}
std::vector<double> ycoords(ny);
for (size_t i = 0; i < ny; ++i) {
ycoords[i] = drawer.minPt().y + i * params.gridResolution;
}
std::unique_ptr<double[]> grid(new double[nx * ny]);
// populate the grid from the gaussians:
for (size_t ix = 0; ix < nx; ++ix) {
auto px = drawer.minPt().x + ix * params.gridResolution;
for (size_t iy = 0; iy < ny; ++iy) {
auto py = drawer.minPt().y + iy * params.gridResolution;
Point2D pt(px, py);
double accum = 0.0;
for (size_t ig = 0; ig < locs.size(); ++ig) {
auto d2 = (pt - locs[ig]).lengthSq();
auto contrib = weights[ig] / widths[ig] *
exp(-0.5 * d2 / (widths[ig] * widths[ig]));
accum += contrib;
}
grid[ix * ny + iy] = accum / (2 * M_PI);
}
}
// and render it:
ContourParams paramsCopy = params;
paramsCopy.setScale = false; // if scaling was needed, we did it already
contourAndDrawGrid(drawer, grid.get(), xcoords, ycoords, nContours, levels,
paramsCopy);
};
// ****************************************************************************
void drawMolACS1996(MolDraw2D &drawer, const ROMol &mol,
const std::string &legend,
const std::vector<int> *highlight_atoms,
const std::vector<int> *highlight_bonds,
const std::map<int, DrawColour> *highlight_atom_map,
const std::map<int, DrawColour> *highlight_bond_map,
const std::map<int, double> *highlight_radii, int confId) {
if (drawer.width() != -1 || drawer.height() != -1) {
BOOST_LOG(rdWarningLog)
<< "ACS drawing mode works best with a flexiCanvas i.e. a drawer"
<< " created with width and height of -1. The scale will be fixed,"
<< " and that may not look great with a pre-determined size."
<< std::endl;
}
auto setAndGo = [&](const ROMol &theMol) -> void {
auto meanBondLen = MolDraw2DUtils::meanBondLength(theMol, confId);
setACS1996Options(drawer.drawOptions(), meanBondLen);
drawer.drawMolecule(theMol, legend, highlight_atoms, highlight_bonds,
highlight_atom_map, highlight_bond_map, highlight_radii,
confId);
};
if (!mol.getNumConformers()) {
// compute 2D coordinates in a standard orientation. This needs to be
// done on a copy because mol is const.
const bool canonOrient = true;
RWMol cpy(mol);
RDDepict::compute2DCoords(cpy, nullptr, canonOrient);
setAndGo(cpy);
} else {
setAndGo(mol);
}
}
// ****************************************************************************
void setACS1996Options(MolDrawOptions &opts, double meanBondLen) {
if (meanBondLen <= 0.0) {
throw ValueErrorException(
"ACS1996Options requires mean bond length > 0.0.");
}
opts.bondLineWidth = 0.6;
opts.scaleBondWidth = false;
// the guideline is for a bond length of 14.4px, and we set things up
// in pixels per Angstrom.
opts.scalingFactor = 14.4 / meanBondLen;
// setting the fixedBondLength means the drawing won't be scaled
// up in a drawer of defined size, so the bond length won't exceed
// 14.4 pixels.
opts.fixedBondLength = 14.4 / meanBondLen;
// offset for multiple bonds is 18% of the bond length.
opts.multipleBondOffset = 0.18;
opts.highlightBondWidthMultiplier = 32;
setMonochromeMode(opts, DrawColour(0.0, 0.0, 0.0), DrawColour(1.0, 1.0, 1.0));
opts.fixedFontSize = 10;
opts.additionalAtomLabelPadding = 0.066;
// The guidelines say Arial font, which is not a free font. A close
// approximation is FreeSans, but that is under GPL v3.0, so can't be
// embedded. Use it if it's there, but fall back on the Roboto font
// which uses an Apache 2.0 license and is also fairly close to Arial.
// It is up to the user to put the FreeSans.ttf in the right place.
// If the user has already specified a fontFile, assume they know
// what they're doing and use it.
if (opts.fontFile.empty()) {
const char *rdbase = getenv("RDBASE");
bool have_free_sans = false;
if (rdbase) {
opts.fontFile = std::string(rdbase) + "/Data/Fonts/FreeSans.ttf";
struct stat buffer;
have_free_sans = (stat(opts.fontFile.c_str(), &buffer) == 0);
}
if (!rdbase || !have_free_sans) {
opts.fontFile = "BuiltinRobotoRegular";
}
}
}
// ****************************************************************************
double meanBondLength(const ROMol &mol, int confId) {
double bondLen = 0.0;
if (mol.getNumBonds()) {
auto conf = mol.getConformer(confId);
for (auto bond : mol.bonds()) {
bondLen += MolTransforms::getBondLength(conf, bond->getBeginAtomIdx(),
bond->getEndAtomIdx());
}
bondLen /= mol.getNumBonds();
}
return bondLen;
}
} // namespace MolDraw2DUtils
} // namespace RDKit
Reference in New Issue View Git Blame Copy Permalink