mirror of
https://github.com/rdkit/rdkit.git
synced 2026-06-03 21:44:30 +08:00
d8f4afb558
Adds a new function MolDraw2D_detail::getSGroupDataLabels() that returns the text and molecule-coordinate positions of DAT SGroup labels, using the same placement logic as the drawing code. This allows external renderers to display SGroup labels consistently with RDKit's placement. Refactors DrawMol::extractSGroupData() to call getSGroupDataLabels() internally, eliminating the duplicate FIELDDISP parsing and position computation logic. Closes #7829 Co-authored-by: Claude Sonnet 4.6 <noreply@anthropic.com>
529 lines
17 KiB
C++
529 lines
17 KiB
C++
//
|
|
// Copyright (C) 2015-2020 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.
|
|
//
|
|
|
|
#include <GraphMol/MolDraw2D/MolDraw2DDetails.h>
|
|
#include <GraphMol/MolDraw2D/StringRect.h>
|
|
#include <GraphMol/Chirality.h>
|
|
#include <GraphMol/FileParsers/FileParserUtils.h>
|
|
#include <GraphMol/SubstanceGroup.h>
|
|
#include <Geometry/Transform2D.h>
|
|
#include <RDGeneral/types.h>
|
|
|
|
#include <cmath>
|
|
#ifndef M_PI
|
|
#define M_PI 3.14159265358979323846
|
|
#endif
|
|
|
|
// ****************************************************************************
|
|
|
|
namespace RDKit {
|
|
namespace MolDraw2D_detail {
|
|
// implementation from $RDBASE/rdkit/sping/pid.py
|
|
void arcPoints(const Point2D &cds1, const Point2D &cds2,
|
|
std::vector<Point2D> &res, float startAng, float extent) {
|
|
// Note: this implementation is simple and not particularly efficient.
|
|
float xScale = (cds2.x - cds1.x) / 2.0;
|
|
float yScale = (cds2.y - cds1.y) / 2.0;
|
|
if (xScale < 0) {
|
|
xScale *= -1;
|
|
}
|
|
if (yScale < 0) {
|
|
yScale *= -1;
|
|
}
|
|
|
|
float x = std::min(cds1.x, cds2.x) + xScale;
|
|
float y = std::min(cds1.y, cds2.y) + yScale;
|
|
|
|
int steps = std::max(static_cast<int>(extent * 2), 5);
|
|
float step = M_PI * extent / (180 * steps);
|
|
float angle = M_PI * startAng / 180;
|
|
for (int i = 0; i <= steps; ++i) {
|
|
Point2D point(x + xScale * cos(angle), y - yScale * sin(angle));
|
|
res.emplace_back(point);
|
|
angle += step;
|
|
}
|
|
}
|
|
|
|
namespace {
|
|
// note, this is approximate since we're just using it for drawing
|
|
bool lineSegmentsIntersect(const Point2D &s1, const Point2D &s2,
|
|
const Point2D &s3, const Point2D &s4) {
|
|
auto d1x = (s1.x - s2.x);
|
|
auto d1y = (s1.y - s2.y);
|
|
auto d2x = (s3.x - s4.x);
|
|
auto d2y = (s3.y - s4.y);
|
|
|
|
if (fabs(d1x) < 1e-4) {
|
|
// fudge factor, since this isn't super critical
|
|
d1x = 1e-4;
|
|
}
|
|
if (fabs(d2x) < 1e-4) {
|
|
// fudge factor, since this isn't super critical
|
|
d2x = 1e-4;
|
|
}
|
|
|
|
auto m1 = d1y / d1x;
|
|
auto m2 = d2y / d2x;
|
|
if (m1 == m2 || m1 == -m2) {
|
|
// parallel
|
|
return false;
|
|
}
|
|
auto b1 = (s1.x * s2.y - s2.x * s1.y) / d1x;
|
|
auto b2 = (s3.x * s4.y - s4.x * s3.y) / d2x;
|
|
|
|
auto intersectX = (b2 - b1) / (m1 - m2);
|
|
return ((intersectX < s1.x) ^ (intersectX < s2.x)) &&
|
|
((intersectX < s3.x) ^ (intersectX < s4.x));
|
|
}
|
|
} // namespace
|
|
|
|
std::vector<Point2D> getBracketPoints(
|
|
const Point2D &p1, const Point2D &p2, const Point2D &refPt,
|
|
const std::vector<std::pair<Point2D, Point2D>> &bondSegments,
|
|
double bracketFrac) {
|
|
std::vector<Point2D> res;
|
|
auto v = p2 - p1;
|
|
Point2D bracketDir{v.y, -v.x};
|
|
bracketDir *= bracketFrac;
|
|
|
|
// we'll default to use the refPt
|
|
auto refVect = p2 - refPt;
|
|
// but check if we intersect any of the bonds:
|
|
for (const auto &seg : bondSegments) {
|
|
if (lineSegmentsIntersect(p1, p2, seg.first, seg.second)) {
|
|
refVect = p2 - seg.first;
|
|
}
|
|
}
|
|
if (bracketDir.dotProduct(refVect) > 0) {
|
|
bracketDir *= -1;
|
|
}
|
|
auto p0 = p1 + bracketDir;
|
|
auto p3 = p2 + bracketDir;
|
|
return {p0, p1, p2, p3};
|
|
}
|
|
|
|
// there are a several empirically determined constants here.
|
|
std::vector<Point2D> handdrawnLine(Point2D cds1, Point2D cds2, double scale,
|
|
bool shiftBegin, bool shiftEnd,
|
|
unsigned nSteps, double deviation,
|
|
double endShift) {
|
|
// std::cout << " " << scale << " " << endShift / scale << std::endl;
|
|
while (endShift / scale > 0.02) {
|
|
endShift *= 0.75;
|
|
}
|
|
if (shiftBegin) {
|
|
cds1.x += (std::rand() % 10 >= 5 ? endShift : -endShift) / scale;
|
|
cds1.y += (std::rand() % 10 >= 5 ? endShift : -endShift) / scale;
|
|
}
|
|
if (shiftEnd) {
|
|
cds2.x += (std::rand() % 10 >= 5 ? endShift : -endShift) / scale;
|
|
cds2.y += (std::rand() % 10 >= 5 ? endShift : -endShift) / scale;
|
|
}
|
|
|
|
Point2D step = (cds2 - cds1) / nSteps;
|
|
// make sure we aren't adding loads of wiggles to short lines
|
|
while (step.length() < 0.2 && nSteps > 2) {
|
|
--nSteps;
|
|
step = (cds2 - cds1) / nSteps;
|
|
}
|
|
// make sure the wiggles aren't too big
|
|
while (deviation / step.length() > 0.15 || deviation * scale > 0.70) {
|
|
deviation *= 0.75;
|
|
}
|
|
Point2D perp{step.y, -step.x};
|
|
perp.normalize();
|
|
std::vector<Point2D> pts;
|
|
pts.push_back(cds1);
|
|
for (unsigned int i = 1; i < nSteps; ++i) {
|
|
auto tgt = cds1 + step * i;
|
|
tgt += perp * deviation * (std::rand() % 20 - 10) / 10.0;
|
|
pts.push_back(tgt);
|
|
}
|
|
pts.push_back(cds2);
|
|
return pts;
|
|
}
|
|
|
|
// ****************************************************************************
|
|
bool doesLineIntersect(const StringRect &rect, const Point2D &end1,
|
|
const Point2D &end2, double padding) {
|
|
Point2D tl, tr, bl, br;
|
|
rect.calcCorners(tl, tr, br, bl, padding);
|
|
if (doLinesIntersect(end2, end1, tl, tr, nullptr)) {
|
|
return true;
|
|
}
|
|
if (doLinesIntersect(end2, end1, tr, br, nullptr)) {
|
|
return true;
|
|
}
|
|
if (doLinesIntersect(end2, end1, br, bl, nullptr)) {
|
|
return true;
|
|
}
|
|
if (doLinesIntersect(end2, end1, bl, tl, nullptr)) {
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
|
|
// ****************************************************************************
|
|
bool doesTriangleIntersect(const StringRect &rect, const Point2D &pt1,
|
|
const Point2D &pt2, const Point2D &pt3,
|
|
double padding) {
|
|
// the quick test is for any of the triangle points inside the rectangle.
|
|
if (rect.isPointInside(pt1, padding) || rect.isPointInside(pt2, padding) ||
|
|
rect.isPointInside(pt3, padding)) {
|
|
return true;
|
|
}
|
|
// But if the rectangle is inside the triangle, that's not enough of a test.
|
|
Point2D tl, tr, br, bl;
|
|
rect.calcCorners(tl, tr, br, bl, padding);
|
|
if (isPointInTriangle(tl, pt1, pt2, pt3) ||
|
|
isPointInTriangle(tr, pt1, pt2, pt3) ||
|
|
isPointInTriangle(br, pt1, pt2, pt3) ||
|
|
isPointInTriangle(bl, pt1, pt2, pt3)) {
|
|
return true;
|
|
}
|
|
// and finally all the points in the rectangle can be outside the triangle,
|
|
// but the sides can cross it. And vice versa. So see if any of the sides
|
|
// intersect.
|
|
if (doLinesIntersect(tl, tr, pt1, pt2, nullptr) ||
|
|
doLinesIntersect(tl, tr, pt2, pt3, nullptr) ||
|
|
doLinesIntersect(tl, tr, pt3, pt1, nullptr) ||
|
|
doLinesIntersect(tr, br, pt1, pt2, nullptr) ||
|
|
doLinesIntersect(tr, br, pt2, pt3, nullptr) ||
|
|
doLinesIntersect(tr, br, pt3, pt1, nullptr) ||
|
|
doLinesIntersect(br, bl, pt1, pt2, nullptr) ||
|
|
doLinesIntersect(br, bl, pt2, pt3, nullptr) ||
|
|
doLinesIntersect(br, bl, pt3, pt1, nullptr) ||
|
|
doLinesIntersect(bl, tl, pt1, pt2, nullptr) ||
|
|
doLinesIntersect(bl, tl, pt2, pt3, nullptr) ||
|
|
doLinesIntersect(bl, tl, pt3, pt1, nullptr)) {
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
|
|
// ****************************************************************************
|
|
bool doesLineIntersectEllipse(const Point2D ¢re, double xradius,
|
|
double yradius, double padding,
|
|
const Point2D &end1, const Point2D &end2) {
|
|
// using
|
|
// https://math.stackexchange.com/questions/76457/check-if-a-point-is-within-an-ellipse
|
|
// to see if either end1 or end2 are inside the ellipse.
|
|
double xr2 = (xradius + padding) * (xradius + padding);
|
|
double yr2 = (yradius + padding) * (yradius + padding);
|
|
double xdisc = (end1.x - centre.x) * (end1.x - centre.x) / xr2;
|
|
double ydisc = (end1.y - centre.y) * (end1.y - centre.y) / yr2;
|
|
if (xdisc + ydisc <= 1.0) {
|
|
return true;
|
|
}
|
|
xdisc = (end2.x - centre.x) * (end2.x - centre.x) / xr2;
|
|
ydisc = (end2.y - centre.y) * (end2.y - centre.y) / yr2;
|
|
return xdisc + ydisc <= 1.0;
|
|
}
|
|
|
|
// ****************************************************************************
|
|
bool doesLineIntersectArc(const Point2D ¢re, double xradius, double yradius,
|
|
double start_ang, double stop_ang, double padding,
|
|
const Point2D &end1, const Point2D &end2) {
|
|
double xr = xradius + padding;
|
|
double yr = yradius + padding;
|
|
double xr2 = xr * xr;
|
|
double yr2 = yr * yr;
|
|
|
|
auto pointInArc = [&](const Point2D &p) -> bool {
|
|
double xdisc = p.x * p.x / xr2;
|
|
double ydisc = p.y * p.y / yr2;
|
|
// start_ang can be more than stop_ang, if, for example, the arc goes from
|
|
// 315 to 45.
|
|
if (xdisc + ydisc <= 1.0) {
|
|
// end1 is inside the whole ellipse. See if the angle it makes with the
|
|
// x axis lies between start_and and stop_ang.
|
|
double th = atan2(p.x, p.y) * 180.0 / M_PI;
|
|
if (th < 0.0) {
|
|
th += 360.0;
|
|
}
|
|
if (start_ang < stop_ang) {
|
|
// it's pretty straightforward
|
|
if (th >= start_ang && th <= stop_ang) {
|
|
return true;
|
|
}
|
|
} else {
|
|
// the arc crosses 0.
|
|
if (th >= start_ang && th <= 360.0) {
|
|
return true;
|
|
}
|
|
if (th >= 0.0 && th <= stop_ang) {
|
|
return true;
|
|
}
|
|
}
|
|
}
|
|
return false;
|
|
};
|
|
|
|
Point2D p1 = end1 - centre;
|
|
if (pointInArc(p1)) {
|
|
return true;
|
|
}
|
|
Point2D p2 = end2 - centre;
|
|
return pointInArc(p2);
|
|
}
|
|
|
|
// ****************************************************************************
|
|
bool doLinesIntersect(const Point2D &l1s, const Point2D &l1f,
|
|
const Point2D &l2s, const Point2D &l2f, Point2D *ip) {
|
|
// using spell from answer 2 of
|
|
// https://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect
|
|
double s1_x = l1f.x - l1s.x;
|
|
double s1_y = l1f.y - l1s.y;
|
|
double s2_x = l2f.x - l2s.x;
|
|
double s2_y = l2f.y - l2s.y;
|
|
|
|
double d = (-s2_x * s1_y + s1_x * s2_y);
|
|
if (d == 0.0) {
|
|
// parallel lines.
|
|
return false;
|
|
}
|
|
double s, t;
|
|
s = (-s1_y * (l1s.x - l2s.x) + s1_x * (l1s.y - l2s.y)) / d;
|
|
t = (s2_x * (l1s.y - l2s.y) - s2_y * (l1s.x - l2s.x)) / d;
|
|
|
|
if (s >= 0 && s <= 1 && t >= 0 && t <= 1) {
|
|
if (ip) {
|
|
ip->x = l1s.x + t * s1_x;
|
|
ip->y = l1s.y + t * s1_y;
|
|
}
|
|
return true;
|
|
}
|
|
|
|
return false;
|
|
}
|
|
|
|
// ****************************************************************************
|
|
bool isPointInTriangle(const Point2D &pt, const Point2D &t1, const Point2D &t2,
|
|
const Point2D &t3) {
|
|
double d = ((t2.y - t3.y) * (t1.x - t3.x) + (t3.x - t2.x) * (t1.y - t3.y));
|
|
double a =
|
|
((t2.y - t3.y) * (pt.x - t3.x) + (t3.x - t2.x) * (pt.y - t3.y)) / d;
|
|
double b =
|
|
((t3.y - t1.y) * (pt.x - t3.x) + (t1.x - t3.x) * (pt.y - t3.y)) / d;
|
|
double c = 1 - a - b;
|
|
return 0 <= a && a <= 1 && 0 <= b && b <= 1 && 0 <= c && c <= 1;
|
|
}
|
|
|
|
std::vector<std::tuple<Point2D, Point2D, Point2D, Point2D>> getWavyLineSegments(
|
|
const Point2D &p1, const Point2D &p2, unsigned int nSegments,
|
|
double vertOffset) {
|
|
std::vector<std::tuple<Point2D, Point2D, Point2D, Point2D>> res;
|
|
|
|
PRECONDITION(nSegments > 1, "too few segments");
|
|
|
|
if (nSegments % 2) {
|
|
++nSegments; // we're going to assume an even number of segments
|
|
}
|
|
|
|
Point2D delta = (p2 - p1);
|
|
Point2D perp(delta.y, -delta.x);
|
|
perp.normalize();
|
|
perp *= vertOffset;
|
|
delta /= nSegments;
|
|
|
|
for (unsigned int i = 0; i < nSegments; ++i) {
|
|
Point2D startpt = p1 + delta * i;
|
|
Point2D segpt = startpt + delta;
|
|
Point2D cpt1 = startpt + perp * (i % 2 ? -1 : 1);
|
|
Point2D cpt2 = segpt + perp * (i % 2 ? -1 : 1);
|
|
res.emplace_back(startpt, cpt1, cpt2, segpt);
|
|
}
|
|
return res;
|
|
}
|
|
|
|
RDKIT_MOLDRAW2D_EXPORT void calcArrowHead(Point2D &arrowEnd, Point2D &arrow1,
|
|
Point2D &arrow2,
|
|
const Point2D &arrowBegin,
|
|
double frac, double lineWidth,
|
|
double angle) {
|
|
if (angle < 1.0e-6) {
|
|
return;
|
|
}
|
|
// Allow for the mitre.
|
|
double adjuster = 0.5 * lineWidth / sin(angle);
|
|
double len = (arrowBegin - arrowEnd).length();
|
|
if (len > 1.0e-6) {
|
|
double adjLen = len - adjuster;
|
|
arrowEnd.x = arrowBegin.x + (arrowEnd.x - arrowBegin.x) * adjLen / len;
|
|
arrowEnd.y = arrowBegin.y + (arrowEnd.y - arrowBegin.y) * adjLen / len;
|
|
}
|
|
|
|
auto delta = arrowBegin - arrowEnd;
|
|
double cos_angle = std::cos(angle), sin_angle = std::sin(angle);
|
|
// To have the arrowhead a consistent fraction of the line length, we need
|
|
// the hypotenuse.
|
|
frac /= cos_angle;
|
|
arrow1 = arrowEnd;
|
|
arrow1.x += frac * (delta.x * cos_angle + delta.y * sin_angle);
|
|
arrow1.y += frac * (delta.y * cos_angle - delta.x * sin_angle);
|
|
|
|
arrow2 = arrowEnd;
|
|
arrow2.x += frac * (delta.x * cos_angle - delta.y * sin_angle);
|
|
arrow2.y += frac * (delta.y * cos_angle + delta.x * sin_angle);
|
|
}
|
|
|
|
// ****************************************************************************
|
|
void adjustLineEndForEllipse(const Point2D ¢re, double xradius,
|
|
double yradius, Point2D p1, Point2D &p2) {
|
|
// move everything so the ellipse is centred on the origin.
|
|
p1 -= centre;
|
|
p2 -= centre;
|
|
double a2 = xradius * xradius;
|
|
double b2 = yradius * yradius;
|
|
double A =
|
|
(p2.x - p1.x) * (p2.x - p1.x) / a2 + (p2.y - p1.y) * (p2.y - p1.y) / b2;
|
|
double B = 2.0 * p1.x * (p2.x - p1.x) / a2 + 2.0 * p1.y * (p2.y - p1.y) / b2;
|
|
double C = p1.x * p1.x / a2 + p1.y * p1.y / b2 - 1.0;
|
|
|
|
auto t_to_point = [&](double t) -> Point2D {
|
|
Point2D ret_val;
|
|
ret_val.x = p1.x + (p2.x - p1.x) * t + centre.x;
|
|
ret_val.y = p1.y + (p2.y - p1.y) * t + centre.y;
|
|
return ret_val;
|
|
};
|
|
|
|
double disc = B * B - 4.0 * A * C;
|
|
if (disc < 0.0) {
|
|
// no solutions, leave things as they are. Bit crap, though.
|
|
p2 += centre;
|
|
return;
|
|
} else if (fabs(disc) < 1.0e-6) {
|
|
// 1 solution
|
|
double t = -B / (2.0 * A);
|
|
p2 = t_to_point(t);
|
|
} else {
|
|
// 2 solutions - take the one nearest p1.
|
|
double disc_rt = sqrt(disc);
|
|
double t1 = (-B + disc_rt) / (2.0 * A);
|
|
double t2 = (-B - disc_rt) / (2.0 * A);
|
|
double t;
|
|
// prefer the t between 0 and 1, as that must be between the original
|
|
// points. If both are, prefer the lower, as that will be nearest p1,
|
|
// so on the bit of the ellipse the line comes to first.
|
|
bool t1_ok = (t1 >= 0.0 && t1 <= 1.0);
|
|
bool t2_ok = (t2 >= 0.0 && t2 <= 1.0);
|
|
if (t1_ok && !t2_ok) {
|
|
t = t1;
|
|
} else if (t2_ok && !t1_ok) {
|
|
t = t2;
|
|
} else if (t1_ok && t2_ok) {
|
|
t = std::min(t1, t2);
|
|
} else {
|
|
// the intersections are both outside the line between p1 and p2
|
|
// so don't do anything.
|
|
p2 += centre;
|
|
return;
|
|
}
|
|
p2 = t_to_point(t);
|
|
}
|
|
}
|
|
|
|
// ****************************************************************************
|
|
std::vector<SGroupDataLabel> getSGroupDataLabels(const ROMol &mol,
|
|
double rotate) {
|
|
std::vector<SGroupDataLabel> result;
|
|
if (!mol.getNumConformers()) {
|
|
return result;
|
|
}
|
|
const auto &sgs = getSubstanceGroups(mol);
|
|
if (sgs.empty()) {
|
|
return result;
|
|
}
|
|
|
|
double rot = rotate * M_PI / 180.0;
|
|
RDGeom::Transform2D tform;
|
|
tform.SetTransform(Point2D(0.0, 0.0), rot);
|
|
|
|
const auto &conf = mol.getConformer();
|
|
for (const auto &sg : sgs) {
|
|
std::string typ;
|
|
if (!sg.getPropIfPresent("TYPE", typ) || typ != "DAT") {
|
|
continue;
|
|
}
|
|
|
|
std::string text;
|
|
if (sg.hasProp("DATAFIELDS")) {
|
|
STR_VECT dfs = sg.getProp<STR_VECT>("DATAFIELDS");
|
|
for (const auto &df : dfs) {
|
|
text += df + "|";
|
|
}
|
|
text.pop_back();
|
|
}
|
|
if (text.empty()) {
|
|
continue;
|
|
}
|
|
|
|
int atomIdx = -1;
|
|
if (!sg.getAtoms().empty()) {
|
|
atomIdx = sg.getAtoms()[0];
|
|
}
|
|
|
|
bool located = false;
|
|
std::string fieldDisp;
|
|
Point2D pos(0.0, 0.0);
|
|
if (sg.getPropIfPresent("FIELDDISP", fieldDisp)) {
|
|
double xp =
|
|
FileParserUtils::stripSpacesAndCast<double>(fieldDisp.substr(0, 10));
|
|
double yp =
|
|
FileParserUtils::stripSpacesAndCast<double>(fieldDisp.substr(10, 10));
|
|
// we always invert y for the molecule coords
|
|
pos = Point2D{xp, -yp};
|
|
|
|
if (fieldDisp[25] == 'R') {
|
|
if (atomIdx < 0) {
|
|
// no atom to anchor relative position to, skip
|
|
continue;
|
|
} else if (fabs(xp) > 1e-3 || fabs(yp) > 1e-3) {
|
|
// opposite sign for y
|
|
pos.x += conf.getAtomPos(atomIdx).x;
|
|
pos.y -= conf.getAtomPos(atomIdx).y;
|
|
located = true;
|
|
}
|
|
} else {
|
|
// Absolute position - check for centroid offset set by drawing pipeline
|
|
if (mol.hasProp("_centroidx")) {
|
|
Point2D centroid;
|
|
mol.getProp("_centroidx", centroid.x);
|
|
mol.getProp("_centroidy", centroid.y);
|
|
// opposite sign for y
|
|
pos.x += centroid.x;
|
|
pos.y -= centroid.y;
|
|
}
|
|
located = true;
|
|
}
|
|
tform.TransformPoint(pos);
|
|
}
|
|
|
|
if (!located) {
|
|
if (atomIdx >= 0) {
|
|
const auto &p = conf.getAtomPos(atomIdx);
|
|
pos = Point2D{p.x, p.y};
|
|
} else {
|
|
BOOST_LOG(rdWarningLog)
|
|
<< "FIELDDISP info not found for DAT SGroup which isn't "
|
|
"associated with an atom. SGroup will not be included."
|
|
<< std::endl;
|
|
continue;
|
|
}
|
|
}
|
|
|
|
result.push_back({text, pos, located, atomIdx});
|
|
}
|
|
return result;
|
|
}
|
|
|
|
} // namespace MolDraw2D_detail
|
|
} // namespace RDKit
|