#2412 - User can't correctly save (or make a layout) to RDF/RXN reaction several products or with separate positioned molecules

core/indigo-core/common/math/algebra.h

@@ -23,6 +23,7 @@
 #include <cmath>
 #include <iostream>
 #include <limits>
+#include <vector>
 
 #include "base_c/defs.h"
 #include "base_cpp/exception.h"
@@ -279,6 +280,7 @@ namespace indigo
         DLLEXPORT void rotateL(float si, float co);
         DLLEXPORT void rotateL(Vec2f vec);
         DLLEXPORT void rotateAroundSegmentEnd(const Vec2f& a, const Vec2f& b, float angle);
+        DLLEXPORT float relativeCross(const Vec2f& a, const Vec2f& b);
 
         DLLEXPORT static float distSqr(const Vec2f& a, const Vec2f& b);
         DLLEXPORT static float dist(const Vec2f& a, const Vec2f& b);
@@ -822,5 +824,234 @@ namespace indigo
         float _d;
     };
 
+    inline bool isPointInPolygon(const Vec2f& p, const std::vector<Vec2f>& poly)
+    {
+        bool in = false;
+        for (size_t i = 0, n = poly.size(); i < n; ++i)
+        {
+            size_t j = (i + 1) % n;
+            bool intersect =
+                ((poly[i].y > p.y) != (poly[j].y > p.y)) && (p.x < (poly[j].x - poly[i].x) * (p.y - poly[i].y) / (poly[j].y - poly[i].y) + poly[i].x);
+            if (intersect)
+                in = !in;
+        }
+        return in;
+    }
+
+    inline std::vector<Vec2f> getPointsInsidePolygon(const std::vector<Vec2f>& polygon1, const std::vector<Vec2f>& polygon2)
+    {
+        std::vector<Vec2f> result;
+        result.reserve(polygon1.size());
+        std::copy_if(polygon1.begin(), polygon1.end(), std::back_inserter(result), [&](auto& p) { return isPointInPolygon(p, polygon2); });
+        return result;
+    }
+
+    inline bool isPointInConvexPolygon(const Vec2f& p, const std::vector<Vec2f>& poly)
+    {
+        auto cross = [](const Vec2f& a, const Vec2f& b, const Vec2f& c) { return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x); };
+        bool sign = cross(poly.back(), poly[0], p) < 0;
+        for (size_t i = 0, n = poly.size(); i < n; ++i)
+            if ((cross(poly[i], poly[(i + 1) % n], p) < 0) != sign)
+                return false;
+        return true;
+    }
+
+    inline std::vector<Vec2f> getPointsInsideConvexPolygon(const std::vector<Vec2f>& polygon1, const std::vector<Vec2f>& polygon2)
+    {
+        std::vector<Vec2f> result;
+        result.reserve(polygon1.size());
+        std::copy_if(polygon1.begin(), polygon1.end(), std::back_inserter(result), [&](auto& p) { return isPointInConvexPolygon(p, polygon2); });
+        return result;
+    }
+
+    inline float convexPolygonArea(const std::vector<Vec2f>& poly)
+    {
+        float area = 0.0f;
+        size_t n = poly.size();
+        for (size_t i = 0; i < n; ++i)
+        {
+            const Vec2f& p1 = poly[i];
+            const Vec2f& p2 = poly[(i + 1) % n];
+            area += Vec2f::cross(p1, p2);
+        }
+        return std::abs(area) * 0.5f;
+    }
+
+    inline bool isInside(const Vec2f& edgeStart, const Vec2f& edgeEnd, const Vec2f& point)
+    {
+        return (edgeEnd.x - edgeStart.x) * (point.y - edgeStart.y) - (edgeEnd.y - edgeStart.y) * (point.x - edgeStart.x) <= 0;
+    }
+
+    inline bool convexPolygonsIntersect(const std::vector<Vec2f>& poly1, const std::vector<Vec2f>& poly2)
+    {
+        auto project = [](const std::vector<Vec2f>& poly, const Vec2f& axis) {
+            auto [min, max] = std::minmax_element(poly.begin(), poly.end(), [&](const Vec2f& a, const Vec2f& b) { return a * axis < b * axis; });
+            return std::pair{*min * axis, *max * axis};
+        };
+
+        auto overlap = [](const auto& range1, const auto& range2) { return !(range1.first > range2.second || range2.first > range1.second); };
+
+        auto getNormals = [](const std::vector<Vec2f>& poly) {
+            std::vector<Vec2f> normals;
+            for (size_t i = 0; i < poly.size(); ++i)
+            {
+                Vec2f edge = poly[(i + 1) % poly.size()] - poly[i];
+                normals.emplace_back(-edge.y, edge.x);
+            }
+            return normals;
+        };
+
+        auto normals1 = getNormals(poly1), normals2 = getNormals(poly2);
+        for (const auto& normal : normals1)
+            if (!overlap(project(poly1, normal), project(poly2, normal)))
+                return false;
+        for (const auto& normal : normals2)
+            if (!overlap(project(poly1, normal), project(poly2, normal)))
+                return false;
+        return true;
+    }
+
+    inline Vec2f computeIntersection(const Vec2f& s, const Vec2f& e, const Vec2f& edgeStart, const Vec2f& edgeEnd)
+    {
+        float dx1 = e.x - s.x, dy1 = e.y - s.y, dx2 = edgeEnd.x - edgeStart.x, dy2 = edgeEnd.y - edgeStart.y;
+        float det = dx1 * dy2 - dy1 * dx2;
+        if (std::abs(det) < 1e-6)
+            return s;
+        float t = ((edgeStart.x - s.x) * dy2 - (edgeStart.y - s.y) * dx2) / det;
+        return {s.x + t * dx1, s.y + t * dy1};
+    }
+
+    inline std::vector<Vec2f> convexClip(const std::vector<Vec2f>& subject, const std::vector<Vec2f>& clip)
+    {
+        std::vector<Vec2f> result = subject;
+        for (size_t i = 0; i < clip.size(); ++i)
+        {
+            const Vec2f& edgeStart = clip[i];
+            const Vec2f& edgeEnd = clip[(i + 1) % clip.size()];
+            std::vector<Vec2f> input = std::move(result);
+            result.clear();
+            if (input.empty())
+                break;
+            Vec2f prev = input.back();
+            for (const Vec2f& curr : input)
+            {
+                if (isInside(edgeStart, edgeEnd, curr))
+                {
+                    if (!isInside(edgeStart, edgeEnd, prev))
+                        result.push_back(computeIntersection(prev, curr, edgeStart, edgeEnd));
+                    result.push_back(curr);
+                }
+                else if (isInside(edgeStart, edgeEnd, prev))
+                {
+                    result.push_back(computeIntersection(prev, curr, edgeStart, edgeEnd));
+                }
+                prev = curr;
+            }
+        }
+        return result;
+    }
+
+    inline float computeConvexContainment(const std::vector<Vec2f>& poly1, const std::vector<Vec2f>& poly2)
+    {
+        float area = convexPolygonArea(poly1);
+        if (area == 0.0f)
+            return 0.0f;
+        std::vector<Vec2f> intersection = convexClip(poly1, poly2);
+        if (intersection.empty())
+            return 0.0f;
+        float intersectionArea = convexPolygonArea(intersection);
+        if (std::abs(intersectionArea - area) < 1e-6f)
+            return 1.0f;
+        return intersectionArea / area;
+    }
+
+    inline float distancePointToEdge(const Vec2f& p, const Vec2f& a, const Vec2f& b)
+    {
+        Vec2f ab = b - a;
+        Vec2f ap = p - a;
+        float t = Vec2f::dot(ap, ab) / Vec2f::dot(ab, ab);
+        t = std::max(0.0f, std::min(1.0f, t));
+        Vec2f projection(a.x + ab.x * t, a.y + ab.y * t);
+        Vec2f diff = p - projection;
+        return diff.length();
+    }
+
+    inline float computeConvexDistance(const std::vector<Vec2f>& poly1, const std::vector<Vec2f>& poly2)
+    {
+        float minDist = std::numeric_limits<float>::max();
+        for (const auto& p : poly1)
+        {
+            float dist = std::numeric_limits<float>::max();
+            size_t n = poly2.size();
+            for (size_t i = 0; i < n; ++i)
+            {
+                float d = distancePointToEdge(p, poly2[i], poly2[(i + 1) % n]);
+                dist = std::min(dist, d);
+                if (d < minDist)
+                    break;
+            }
+            minDist = std::min(minDist, dist);
+        }
+        for (const auto& p : poly2)
+        {
+            float dist = std::numeric_limits<float>::max();
+            size_t n = poly1.size();
+            for (size_t i = 0; i < n; ++i)
+            {
+                float d = distancePointToEdge(p, poly1[i], poly1[(i + 1) % n]);
+                dist = std::min(dist, d);
+                if (d < minDist)
+                    break;
+            }
+            minDist = std::min(minDist, dist);
+        }
+        return minDist;
+    }
+
+    inline bool doesVerticalLineIntersectPolygon(float x, const std::vector<Vec2f>& poly)
+    {
+        for (size_t i = 0, n = poly.size(); i < n; ++i)
+        {
+            const Vec2f& a = poly[i];
+            const Vec2f& b = poly[(i + 1) % n];
+            if ((x > std::min(a.x, b.x) && x < std::max(a.x, b.x)))
+                return true;
+        }
+        return false;
+    }
+
+    inline bool doesHorizontalLineIntersectPolygon(float y, const std::vector<Vec2f>& poly)
+    {
+        for (size_t i = 0, n = poly.size(); i < n; ++i)
+        {
+            const Vec2f& a = poly[i];
+            const Vec2f& b = poly[(i + 1) % n];
+            if ((y > std::min(a.y, b.y) && y < std::max(a.y, b.y)))
+                return true;
+        }
+        return false;
+    }
+
+    inline bool doesRayIntersectPolygon(const Vec2f& p1, const Vec2f& p2, const std::vector<Vec2f>& poly)
+    {
+        Vec2f dir = p2 - p1;
+        for (size_t i = 0, n = poly.size(); i < n; ++i)
+        {
+            Vec2f A = poly[i];
+            Vec2f B = poly[(i + 1) % n];
+            Vec2f a = A - p1;
+            Vec2f b = B - p1;
+            Vec2f s = b - a;
+            float den = Vec2f::cross(dir, s);
+            if (std::fabs(den) <= +0.0f)
+                continue;
+            float t = Vec2f::cross(a, s) / den;
+            float u = Vec2f::cross(a, dir) / den;
+            if (t >= 0.f && u >= 0.f && u <= 1.f)
+                return true;
+        }
+        return false;
+    }
+
 } // namespace indigo
 #endif

core/indigo-core/common/math/vec2f.cpp

@@ -189,6 +189,11 @@ float Vec2f::cross(const Vec2f& a, const Vec2f& b)
     return a.x * b.y - a.y * b.x;
 }
 
+float Vec2f::relativeCross(const Vec2f& a, const Vec2f& b)
+{
+    return (a.x - x) * (b.y - y) - (a.y - y) * (b.x - x);
+}
+
 void Vec2f::projectZ(Vec2f& v2, const Vec3f& v3)
 {
     v2.x = v3.x;

core/indigo-core/molecule/base_molecule.h

@@ -589,6 +589,9 @@ namespace indigo
         void getBoundingBox(float font_size, LABEL_MODE label_mode, Vec2f& bottom_left, Vec2f& top_right);
         void getBoundingBox(float font_size, LABEL_MODE label_mode, Rect2f& bbox);
 
+        // calc convex hull
+        std::vector<Vec2f> getConvexHull(const Vec2f& min_box) const;
+
         // aliases
         bool isAlias(int atom_idx) const;
         const char* getAlias(int atom_idx) const;

core/indigo-core/molecule/src/base_molecule.cpp

@@ -4775,6 +4775,40 @@ void BaseMolecule::getBoundingBox(float font_size, LABEL_MODE label_mode, Rect2f
     bbox = Rect2f(a, b);
 }
 
+std::vector<Vec2f> BaseMolecule::getConvexHull(const Vec2f& min_box) const
+{
+    std::vector<Vec2f> vertices;
+    std::transform(_xyz.ptr(), _xyz.ptr() + _xyz.size(), std::back_inserter(vertices), [](const Vec3f& v) -> Vec2f { return Vec2f(v.x, v.y); });
+    if (vertices.size() < 3)
+    {
+        Rect2f bbox;
+        getBoundingBox(bbox, min_box);
+        vertices.clear();
+        vertices.emplace_back(bbox.leftTop());
+        vertices.emplace_back(bbox.rightTop());
+        vertices.emplace_back(bbox.rightBottom());
+        vertices.emplace_back(bbox.leftBottom());
+        return vertices;
+    }
+    std::sort(vertices.begin(), vertices.end());
+    std::vector<Vec2f> hull;
+    for (const auto& p : vertices)
+    {
+        while (hull.size() >= 2 && hull[hull.size() - 2].relativeCross(hull.back(), p) <= 0)
+            hull.pop_back();
+        hull.push_back(p);
+    }
+    size_t lower_size = hull.size();
+    for (auto it = vertices.rbegin(); it != vertices.rend(); ++it)
+    {
+        while (hull.size() > lower_size && hull[hull.size() - 2].relativeCross(hull.back(), *it) <= 0)
+            hull.pop_back();
+        hull.push_back(*it);
+    }
+    hull.pop_back();
+    return hull;
+}
+
 void BaseMolecule::getBoundingBox(Vec2f& a, Vec2f& b) const
 {
     for (int atom_idx = 0; atom_idx < vertexCount(); ++atom_idx)

core/indigo-core/reaction/reaction_multistep_detector.h

@@ -29,6 +29,7 @@
 #include <vector>
 
 #include "base_cpp/exception.h"
+#include "layout/metalayout.h"
 #include "molecule/meta_commons.h"
 
 namespace indigo
@@ -41,13 +42,36 @@ namespace indigo
     class ReactionMultistepDetector
     {
     public:
+        using MOL_DISTANCES = std::vector<std::pair<size_t, float>>;
+        using MOL_DISTANCES_MAP = std::unordered_map<size_t, float>;
+
+        struct MOL_DISTANCES_DESC
+        {
+            MOL_DISTANCES sorted_distances;
+            MOL_DISTANCES_MAP distances_map;
+        };
+
         enum class ReactionType
         {
             ESimpleReaction,
             EMutistepReaction,
             EPathwayReaction
         };
 
+        enum class ZoneType
+        {
+            EPlus,
+            EArrow,
+            EPathWay
+        };
+
+        struct SPECIAL_ZONE_DESC
+        {
+            ZoneType zone_type;
+            std::vector<std::vector<Vec2f>> zone_sections;
+            std::vector<Vec2f> origin_coordinates;
+        };
+
         ReactionMultistepDetector(BaseMolecule& mol);
         ~ReactionMultistepDetector();
         ReactionType detectReaction();
@@ -62,11 +86,25 @@ namespace indigo
         typedef std::vector<FLOAT_INT_PAIR> FLOAT_INT_PAIRS;
         const Vec2f PLUS_BBOX_SHIFT = {0.9f, 0.9f};
         const Vec2f ARROW_BBOX_SHIFT = {0.0f, 0.9f};
+        const float PLUS_DETECTION_DISTANCE = LayoutOptions::DEFAULT_BOND_LENGTH * 2.5;
+        const float ARROW_DETECTION_DISTANCE = LayoutOptions::DEFAULT_BOND_LENGTH * 3;
 
         DECL_ERROR;
 
     private:
         void createSummBlocks();
+        // collect molecules' distances
+        void collectSortedDistances();
+        void createSpecialZones();
+        void addPlusZones(const Vec2f& pos);
+        void addArrowZones(const Vec2f& tail, const Vec2f& head);
+        void addPathwayZones(const Vec2f& head, const Vec2f& sp_beg, const Vec2f& sp_end, const std::vector<Vec2f>& tails);
+        std::map<int, std::unordered_set<int>> findSpecialZones(size_t mol_idx);
+        std::optional<std::pair<int, int>> findMaxSpecialZone(size_t mol_idx);
+        void mergeCloseComponents();
+        std::optional<std::pair<int, int>> isMergeable(size_t mol_idx1, size_t mol_idx2, std::optional<std::pair<int, int>> current_zone);
+        bool checkForOppositeSections(ZoneType zt, const std::unordered_set<int>& sections1, const std::unordered_set<int>& sections2);
+        std::unique_ptr<BaseMolecule> extractComponent(int index);
         void sortSummblocks();
 
         bool mapReactionComponents();
@@ -79,7 +117,9 @@ namespace indigo
         std::vector<ReactionComponent> _reaction_components;
         std::vector<MolSumm> _component_summ_blocks;
         std::list<MolSumm> _component_summ_blocks_list;
-
+        std::vector<std::pair<std::unique_ptr<BaseMolecule>, std::vector<Vec2f>>> _components;
+        std::vector<MOL_DISTANCES_DESC> _mol_distances;
+        std::vector<SPECIAL_ZONE_DESC> _zones;
         int _moleculeCount;
     };