C++ : Merging two convex hulls

Question

I wrote a C++ program that merges two convex hulls in linear time. The algorithm is based on http://cgm.cs.mcgill.ca/~godfried/teaching/cg-projects/97/Plante/CompGeomProject-EPlante/algorithm.html

Feel free to comment anything!

Point.h

#include <cassert>
#include <cmath>
#include <iostream>

constexpr double tolerance = 1e-6;

struct Point2d {
    double x = 0.0;
    double y = 0.0;
    Point2d() = default;
    Point2d(double x, double y) : x {x}, y {y} {}
};

bool isClose(double x, double d) {
    return std::fabs(x - d) < tolerance;
}

std::ostream& operator<<(std::ostream& os, const Point2d& point) {
    os << '{' << point.x << ", " << point.y << '}';
    return os;
}

Point2d operator+(const Point2d& lhs, const Point2d& rhs) {
    return Point2d(lhs.x + rhs.x, lhs.y + rhs.y);
}

Point2d operator-(const Point2d& lhs, const Point2d& rhs) {
    return Point2d(lhs.x - rhs.x, lhs.y - rhs.y);
}

bool operator==(const Point2d& lhs, const Point2d& rhs) {
    return isClose(lhs.x, rhs.x) && isClose(lhs.y, rhs.y);
}

bool operator!=(const Point2d& lhs, const Point2d& rhs) {
    return !(lhs == rhs);
}

bool operator<(const Point2d& lhs, const Point2d& rhs) {
    return rhs.x - lhs.x > tolerance ||
           (isClose(lhs.x, rhs.x) && rhs.y - lhs.y > tolerance);
}

bool operator>=(const Point2d& lhs, const Point2d& rhs) {
    return !(lhs < rhs);
}

bool operator>(const Point2d& lhs, const Point2d& rhs) {
    return lhs.x - rhs.x > tolerance ||
           (isClose(lhs.x, rhs.x) && lhs.y - rhs.y > tolerance);
}

bool operator<=(const Point2d& lhs, const Point2d& rhs) {
    return !(lhs > rhs);
}

double Dot(const Point2d& lhs, const Point2d& rhs) {
    return lhs.x * rhs.x + lhs.y * rhs.y;
}

double Cross(const Point2d& lhs, const Point2d& rhs) {
    return lhs.x * rhs.y - lhs.y * rhs.x;
}

static const Point2d orig {0.0, 0.0};

double dist(const Point2d& lhs, const Point2d& rhs = orig) {
    return std::hypot(lhs.x - rhs.x, lhs.y - rhs.y);
}

double Cos(const Point2d& p1, const Point2d& p2) {
    assert(p1 != orig && p2 != orig);
    auto dot = Dot(p1, p2);
    auto r = dist(p1) * dist(p2);
    return dot / r;
}

double Orientation(const Point2d& p0, const Point2d& p1, const Point2d& p2) {
    auto cross = Cross(p1 - p0, p2 - p0);
    return cross;
}

ConvexHull.h

#include "Point.h"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <iostream>
#include <iterator>
#include <random>
#include <set>
#include <vector>
#include <ranges>

namespace sr = std::ranges;

enum class Order {
    Before,
    EqualSrc,
    Between,
    EqualDst,
    After,
};

Order getOrder(const Point2d& p1, const Point2d& p2, const Point2d& q) {
    assert(p1 != p2);
    if (q == p1) {
        return Order::EqualSrc;
    } else if (q == p2) {
        return Order::EqualDst;
    }
    if (!isClose(p1.x, p2.x)) {
        if (p2.x - p1.x > tolerance) { // p1 - p2
            if (p1.x - q.x > tolerance) { // q - p1 - p2
                return Order::Before;
            } else if (p2.x - q.x > tolerance) { // p1 - q - p2
                return Order::Between;
            } else { // p1 - p2 - q
                return Order::After;
            }
        } else { // p2 - p1
            if (q.x - p1.x > tolerance) { // p2 - p1 - q
                return Order::Before;
            } else if (q.x - p2.x > tolerance) { // p2 - q - p1
                return Order::Between;
            } else { // q - p2 - p1
                return Order::After;
            }
        }
    } else {
        if (p2.y - p1.y > tolerance) { // p1 - p2
            if (p1.y - q.y > tolerance) { // q - p1 - p2
                return Order::Before;
            } else if (p2.y - q.y > tolerance) { // p1 - q - p2
                return Order::Between;
            } else { // p1 - p2 - q
                return Order::After;
            }
        } else { // p2 - p1
            if (q.y - p1.y > tolerance) { // p2 - p1 - q
                return Order::Before;
            } else if (q.y - p2.y > tolerance) { // p2 - q - p1
                return Order::Between;
            } else { // q - p2 - p1
                return Order::After;
            }
        }
    }
}

std::vector<Point2d> ConvexHullAdd(std::vector<Point2d>& CH, const Point2d& p) {
    if (CH.empty()) {
        CH.push_back(p);
        return CH;
    } else if (CH.size() == 1) {
        if (CH[0] != p) {
            CH.push_back(p);
        }
        return CH;
    } else if (CH.size() == 2) {
        const auto& p0 = CH[0];
        const auto& p1 = CH[1];
        auto cross = Cross(p1 - p0, p - p0);
        if (isClose(cross, 0.0)) { // colinear
            auto d1 = dist(p0, p);
            auto d2 = dist(p1, p);
            auto d3 = dist(p0, p1);
            auto d = std::max({d1, d2, d3});
            if (isClose(d, d1)) {
                return std::vector<Point2d>{p0, p};
            } else if (isClose(d, d2)) {
                return std::vector<Point2d>{p1, p};
            } else {
                return std::vector<Point2d>{p0, p1};
            }
        } else if (cross > tolerance) { // left turn
            return std::vector<Point2d>{p0, p, p1};
        } else { // right turn
            return std::vector<Point2d>{p0, p1, p};
        }
    }

    // make clockwise convex hull
    const std::size_t m = CH.size();
    std::size_t i = 0;
    // while right turn
    while (true) {
        if (i >= m) {
            break;
        }
        auto cross = Cross(p - CH[i % m], CH[(i + 1) % m] - CH[i % m]);
        if (isClose(cross, 0.0)) {
            // special case: colinear
            auto ord = getOrder(CH[i % m], CH[(i + 1) % m], p);
            if (ord == Order::Before) {
                CH[i % m] = p;
            } else if (ord == Order::After) {
                CH[(i + 1) % m] = p;
            }
            return CH;
        } else if (cross > tolerance) {
            i++;
        } else {
            break;
        }
    }
    if (i == m) { // interior point
        return CH;
    }
    if (i == 0) {
        // first turn is left, turn back to where left turn begins
        while (true) {
            auto cross = Cross(p - CH[i % m], CH[(i + 1) % m] - CH[i % m]);
            if (isClose(cross, 0.0)) {
                // special case: colinear
                auto ord = getOrder(CH[i % m], CH[(i + 1) % m], p);
                if (ord == Order::Before) {
                    CH[i % m] = p;
                } else if (ord == Order::After) {
                    CH[(i + 1) % m] = p;
                }
                return CH;
            } else if (cross < -tolerance) {
                i = (i + m - 1) % m;
            } else {
                break;
            }
        }
        i = (i + 1) % m;
    }
    std::size_t j = (i + 1) % m;
    // while left turn
    while (true) {
        auto cross = Cross(p - CH[j % m], CH[(j + 1) % m] - CH[j % m]);
        if (isClose(cross, 0.0)) {
            // special case: colinear
            auto ord = getOrder(CH[j % m], CH[(j + 1) % m], p);
            if (ord == Order::Before) {
                CH[j % m] = p;
            } else if (ord == Order::After) {
                CH[(j + 1) % m] = p;
            }
            return CH;
        } else if (cross < -tolerance) {
            j = (j + 1) % m;
        } else {
            break;
        }
    }

    // replace p_(i, j) with p
    if (j < i) {
        CH.erase(CH.begin() + i + 1, CH.end());
        CH.push_back(p);
        CH.erase(CH.begin(), CH.begin() + j);
    } else if (j > i) {
        CH.erase(CH.begin() + i + 1, CH.begin() + j);
        CH.insert(CH.begin() + i + 1, p);
    } else { // cannot happen
        assert(0);
    }
    return CH;
}

using ColinearPair = std::pair<std::optional<Point2d>, std::optional<Point2d>>;

void UpdateColinearPair(ColinearPair& curr_minmax, const Point2d& new_min, const Point2d& new_max) {
    auto& [curr_min, curr_max] = curr_minmax;
    if (curr_min.has_value()) {
        auto ord = getOrder(new_min, new_max, *curr_min);
        if (ord != Order::Before) {
            curr_min = new_min;
        }
    } else {
        curr_min = new_min;
    }
    if (curr_max.has_value()) {
        auto ord = getOrder(new_min, new_max, *curr_max);
        if (ord != Order::After) {
            curr_max = new_max;
        }
    } else {
        curr_max = new_max;
    }
}

void HandleColinearCase(ColinearPair& curr_minmax, const Point2d& p1, const Point2d& p2, const Point2d& q) {
    auto ord = getOrder(p1, p2, q);
    if (ord == Order::Before || ord == Order::EqualSrc) { // (q, p2)
        UpdateColinearPair(curr_minmax, q, p2);
    } else if (ord == Order::Between) { // (p1, p2)
        UpdateColinearPair(curr_minmax, p1, p2);
    } else if (ord == Order::EqualDst || ord == Order::After) { // (p1, q)
        UpdateColinearPair(curr_minmax, p1, q);
    } else { // cannot happen
        assert(0);
    }
}

std::vector<Point2d> MergeTwoConvexHulls(const std::vector<Point2d>& CH1, const std::vector<Point2d>& CH2) {
    if (std::min(CH1.size(), CH2.size()) <= 2) { // base case: just do online update.
        if (CH1.size() < CH2.size()) {
            auto CH = CH2;
            for (const auto& p : CH1) {
                CH = ConvexHullAdd(CH, p);
            }
            return CH;
        } else {
            auto CH = CH1;
            for (const auto& p : CH2) {
                CH = ConvexHullAdd(CH, p);
            }
            return CH;
        }
    }
    const std::size_t m = CH1.size();
    const std::size_t n = CH2.size();
    std::size_t start1 = std::distance(CH1.begin(), sr::min_element(CH1));
    std::size_t start2 = std::distance(CH2.begin(), sr::min_element(CH2));

    std::size_t i = start1;
    std::size_t j = start2;

    // compare angle with vertical axis to decide first edge vs vertex pair
    enum class Edge {
        CH1,
        CH2,
    };

    Edge currEdge = (Cos(Point2d{0.0, 1.0}, CH1[(i + 1) % m] - CH1[i % m])
            >= Cos(Point2d{0.0, 1.0}, CH2[(j + 1) % n] - CH2[j % n]))
                    ? Edge::CH1 : Edge::CH2;

    std::vector<Point2d> CH;

    ColinearPair colinear_minmax;

    while (i < start1 + m || j < start2 + n) {
        const auto& p_i = CH1[i % m];
        const auto& q_j = CH2[j % n];
        const auto& p_i1 = CH1[(i + 1) % m];
        const auto& q_j1 = CH2[(j + 1) % n];
        if (currEdge == Edge::CH1) {
            // compare p_i-p_(i+1) vs q_j
            auto dir = Orientation(p_i, p_i1, q_j);
            if (dir > tolerance) { // q_j is in the left side: add q_j
                CH.push_back(q_j);
            } else if (dir < -tolerance) { // q_j is in the right side: add p_i, p_(i+1)
                CH.push_back(p_i);
                CH.push_back(p_i1);
            } else { // colinear: determine the order
                HandleColinearCase(colinear_minmax, p_i, p_i1, q_j);
            }

            // choose one to advance
            auto curr_edge = p_i1 - p_i;
            const auto& p_i2 = CH1[(i + 2) % m];
            auto next_edge_p = p_i2 - p_i1;
            auto next_edge_q = q_j1 - q_j;

            auto cos_p = Cos(curr_edge, next_edge_p);
            auto cos_q = Cos(curr_edge, next_edge_q);
            if (cos_p - cos_q > 0.0) { // advance p
                i++;
            } else { // advance q
                currEdge = Edge::CH2;
                i++;
            }
            if (isClose(std::max(cos_p, cos_q), 1.0)
                && colinear_minmax.first.has_value()) { // next edge is parallel with curr edge, and in fact colinear
                // do nothing, keep colinear pairs
            } else if (colinear_minmax.first.has_value()) {
                // add this colinear pair to CH
                CH.push_back(*colinear_minmax.first);
                CH.push_back(*colinear_minmax.second);
                colinear_minmax = {{}, {}}; // discard colinear pair
            }
        } else if (currEdge == Edge::CH2) {
            // compare q_j-q_(j+1) vs p_i
            auto dir = Orientation(q_j, q_j1, p_i);
            if (dir > tolerance) { // p_i is in the left side: add p_i
                CH.push_back(p_i);
            } else if (dir < -tolerance) { // p_i is in the right side: add q_j, q_(j+1)
                CH.push_back(q_j);
                CH.push_back(q_j1);
            } else { // colinear: determine the order
                HandleColinearCase(colinear_minmax, q_j, q_j1, p_i);
            }

            // choose one to advance
            auto curr_edge = q_j1 - q_j;
            const auto& q_j2 = CH2[(j + 2) % n];
            auto next_edge_p = p_i1 - p_i;
            auto next_edge_q = q_j2 - q_j1;

            auto cos_p = Cos(curr_edge, next_edge_p);
            auto cos_q = Cos(curr_edge, next_edge_q);
            if (cos_q - cos_p > 0.0) { // advance q
                j++;
            } else { // advance p
                currEdge = Edge::CH1;
                j++;
            }
            if (isClose(std::max(cos_p, cos_q), 1.0)
                && colinear_minmax.first.has_value()) { // next edge is parallel with curr edge, and in fact colinear
                // do nothing, keep colinear pairs
            } else if (colinear_minmax.first.has_value()) {
                // add this colinear pair to CH
                CH.push_back(*colinear_minmax.first);
                CH.push_back(*colinear_minmax.second);
                colinear_minmax = {{}, {}}; // discard colinear pair
            }
        } else { // cannot happen
            assert(0);
        }
    }
    CH.erase(std::unique(CH.begin(), CH.end()), CH.end());
    if (CH.back() == CH.front()) {
        CH.pop_back();
    }
    return CH;
}

Test code:

#include "Point.h"
#include "ConvexHull.h"
#include <iostream>
#include <vector>

int main() {

    // Test case 1: ◇□
    {
        std::vector<Point2d> CH1{{0.0, 0.0},
                                 {2.0, 2.0},
                                 {4.0, 0.0},
                                 {2.0, -2.0}};
        std::vector<Point2d> CH2{{5.0, -1.0},
                                 {5.0, 1.0},
                                 {7.0, 1.0},
                                 {7.0, -1.0}};

        auto CH = MergeTwoConvexHulls(CH1, CH2);
        std::cout << "Convex hull:\n";
        for (const auto& p : CH) {
            std::cout << p << '\n';
        }
        std::cout << '\n';
    }

    // Test case 2: □ □
    {
        std::vector<Point2d> CH1{{0.0, 0.0},
                                 {0.0, 2.0},
                                 {2.0, 2.0},
                                 {2.0, 0.0}};
        std::vector<Point2d> CH2{{4.0, 0.0},
                                 {4.0, 2.0},
                                 {6.0, 2.0},
                                 {6.0, 0.0}};

        auto CH = MergeTwoConvexHulls(CH1, CH2);
        std::cout << "Convex hull:\n";
        for (const auto& p : CH) {
            std::cout << p << '\n';
        }
        std::cout << '\n';
    }

    // Test case 3: □◇ overlapping
    {
        std::vector<Point2d> CH1{{0.0, 0.0},
                                 {0.0, 2.0},
                                 {2.0, 2.0},
                                 {2.0, 0.0}};
        std::vector<Point2d> CH2{{-0.5, 1.0},
                                 {1.0, 2.5},
                                 {2.5, 1.0},
                                 {1.0, -0.5}};

        auto CH = MergeTwoConvexHulls(CH1, CH2);
        std::cout << "Convex hull:\n";
        for (const auto& p : CH) {
            std::cout << p << '\n';
        }
        std::cout << '\n';
    }

    // Test case 4: □ inscribes ◇
    {
        std::vector<Point2d> CH1{{0.0, 0.0},
                                 {0.0, 2.0},
                                 {2.0, 2.0},
                                 {2.0, 0.0}};
        std::vector<Point2d> CH2{{-1.0, 1.0},
                                 {1.0, 3.0},
                                 {3.0, 1.0},
                                 {1.0, -1.0}};

        auto CH = MergeTwoConvexHulls(CH1, CH2);
        std::cout << "Convex hull:\n";
        for (const auto& p : CH) {
            std::cout << p << '\n';
        }
        std::cout << '\n';
    }

    // Test case 5: □□ overlapping
    {
        std::vector<Point2d> CH1{{0.0, 0.0},
                                 {0.0, 2.0},
                                 {2.0, 2.0},
                                 {2.0, 0.0}};
        std::vector<Point2d> CH2{{1.0, 0.0},
                                 {1.0, 2.0},
                                 {3.0, 2.0},
                                 {3.0, 0.0}};

        auto CH = MergeTwoConvexHulls(CH1, CH2);
        std::cout << "Convex hull:\n";
        for (const auto& p : CH) {
            std::cout << p << '\n';
        }
        std::cout << '\n';
    }

    // Test case 6: □ □ meets at corner
    {
        std::vector<Point2d> CH1{{0.0, 0.0},
                                 {0.0, 2.0},
                                 {2.0, 2.0},
                                 {2.0, 0.0}};
        std::vector<Point2d> CH2{{2.0, 2.0},
                                 {2.0, 4.0},
                                 {4.0, 4.0},
                                 {4.0, 2.0}};

        auto CH = MergeTwoConvexHulls(CH1, CH2);
        std::cout << "Convex hull:\n";
        for (const auto& p : CH) {
            std::cout << p << '\n';
        }
        std::cout << '\n';
    }

    // Test case 7: □ contains □, two □ intersects at edges
    {
        std::vector<Point2d> CH1{{0.0, 0.0},
                                 {0.0, 2.0},
                                 {2.0, 2.0},
                                 {2.0, 0.0}};
        std::vector<Point2d> CH2{{0.0, 0.0},
                                 {0.0, 1.0},
                                 {1.0, 1.0},
                                 {1.0, 0.0}};

        auto CH = MergeTwoConvexHulls(CH1, CH2);
        std::cout << "Convex hull:\n";
        for (const auto& p : CH) {
            std::cout << p << '\n';
        }
        std::cout << '\n';
    }

    // Test case 8: □ contains □ with no intersection
    {
        std::vector<Point2d> CH1{{0.0, 0.0},
                                 {0.0, 2.0},
                                 {2.0, 2.0},
                                 {2.0, 0.0}};
        std::vector<Point2d> CH2{{0.5, 0.5},
                                 {0.5, 1.5},
                                 {1.5, 1.5},
                                 {1.5, 0.5}};

        auto CH = MergeTwoConvexHulls(CH1, CH2);
        std::cout << "Convex hull:\n";
        for (const auto& p : CH) {
            std::cout << p << '\n';
        }
        std::cout << '\n';
    }

    // Test case 9: □ and interior point
    {
        std::vector<Point2d> CH1{{0.0, 0.0},
                                 {0.0, 2.0},
                                 {2.0, 2.0},
                                 {2.0, 0.0}};
        std::vector<Point2d> CH2{{1.0, 1.0}};

        auto CH = MergeTwoConvexHulls(CH1, CH2);
        std::cout << "Convex hull:\n";
        for (const auto& p : CH) {
            std::cout << p << '\n';
        }
        std::cout << '\n';
    }

    // Test case 10: □ and exterior point
    {
        std::vector<Point2d> CH1{{0.0, 0.0},
                                 {0.0, 2.0},
                                 {2.0, 2.0},
                                 {2.0, 0.0}};
        std::vector<Point2d> CH2{{3.0, 1.0}};

        auto CH = MergeTwoConvexHulls(CH1, CH2);
        std::cout << "Convex hull:\n";
        for (const auto& p : CH) {
            std::cout << p << '\n';
        }
        std::cout << '\n';
    }

    // Test case 11: □ and boundary point
    {
        std::vector<Point2d> CH1{{0.0, 0.0},
                                 {0.0, 2.0},
                                 {2.0, 2.0},
                                 {2.0, 0.0}};
        std::vector<Point2d> CH2{{2.0, 1.0}};

        auto CH = MergeTwoConvexHulls(CH1, CH2);
        std::cout << "Convex hull:\n";
        for (const auto& p : CH) {
            std::cout << p << '\n';
        }
        std::cout << '\n';
    }

    // Test case 12: □ and exterior point colinear
    {
        std::vector<Point2d> CH1{{0.0, 0.0},
                                 {0.0, 2.0},
                                 {2.0, 2.0},
                                 {2.0, 0.0}};
        std::vector<Point2d> CH2{{2.0, 3.0}};

        auto CH = MergeTwoConvexHulls(CH1, CH2);
        std::cout << "Convex hull:\n";
        for (const auto& p : CH) {
            std::cout << p << '\n';
        }
        std::cout << '\n';
    }

    // Test case 13: □ and exterior point colinear 2
    {
        std::vector<Point2d> CH1{{0.0, 0.0},
                                 {0.0, 2.0},
                                 {2.0, 2.0},
                                 {2.0, 0.0}};
        std::vector<Point2d> CH2{{2.0, -1.0}};

        auto CH = MergeTwoConvexHulls(CH1, CH2);
        std::cout << "Convex hull:\n";
        for (const auto& p : CH) {
            std::cout << p << '\n';
        }
        std::cout << '\n';
    }

    // Test case 14: regular octagon and ▷ that intersects at 4 points
    {
        auto sq2 = std::sqrt(2);
        std::vector<Point2d> CH1{{-2.0, 0.0},
                                 {-sq2, sq2},
                                 {0.0, 2.0},
                                 {sq2, sq2},
                                 {2.0, 0.0},
                                 {sq2, -sq2},
                                 {0.0, -2.0},
                                 {-sq2, -sq2}};
        std::vector<Point2d> CH2{{-2.0, 4.0 - 4.0 * sq2},
                                 {-2.0, -4.0 + 4.0 * sq2},
                                 {2.0, 0.0}};

        auto CH = MergeTwoConvexHulls(CH1, CH2);
        std::cout << "Convex hull:\n";
        for (const auto& p : CH) {
            std::cout << p << '\n';
        }
        std::cout << '\n';
    }

}
````

Answer 1

General impression when starting to read your code is that it's good! You're using constexpr, defining your constructors nicely, etc.

includes

Generally, I'll list standard includes first, and project-local includes at the end of the list.

scope and namespaces

Your point header is defining tolerance as a global variable, which may have nothing to do with other classes or even conflict with other headers. You ought to make it a member of the class, or make everything go inside a namespace.

The various operators are more efficient when declared as "hidden friends" rather than non-members. This makes overload resolution faster and more well behaved when you have complex interactions.

operators

Normally, you would define + in terms of += etc. You defined + only, rather than all related operations.

All the relational operations can be defined, as of C++20, as a single operator<=> function, and some of them can be automatically generated from others rather than you doing the exact same thing explicitly; in particular, operator!= will be automatically generated from operator==.

constexpr

You write

static const Point2d orig {0.0, 0.0};

(which is polluting the global namespace, as indicated in the first point above) why isn't it constexpr? Well, your constructors are not constexpr-friendly. There's no reason why it can't be! Same with many other functions; if the compiler can simplify things at compile-time, why not let it?!

constructors

User indi points out that if you leave off both constructors, you'll automatically get the same result anyway. Thanks to the inline initializers on the fields, the default constructor will apply those, and the newer updated rules for aggregate initializers means you can still use curly-brace syntax like a plain struct. I agree that for a class like this, that's perfectly OK, and would include a comment in lieu of the constructor indicating that it's meant to be used that way.

test code

Your test code is repeated, with only the actual input changing. Your block of work

    auto CH = MergeTwoConvexHulls(CH1, CH2);
    std::cout << "Convex hull:\n";
    for (const auto& p : CH) {
        std::cout << p << '\n';
    }
    std::cout << '\n';

should be a function call. You're even using the same variable names (CH1 and CH2) in every case!

const

std::vector<Point2d> ConvexHullAdd(std::vector<Point2d>& CH, const Point2d& p)

Is the function modifying CH as it works? I don't think it ought to, so define the parameter to be const. I noticed you have repetitive code, not just i%m (expensive division) being repeated but the subscripting of the input vector being repeated: CH[i%m]. Making things const will help ensure that the compiler can eliminate the repeated work automatically.

Answer 2

In addition to JDługosz's excellent answer, I'd like to add these remarks:

Use switch-statements where appropriate

Replace chains of if-else-statements that check the value of an enum type variable with a switch-statement. This usually makes the code more readable, and also has the advantage that the compiler can warn you if you don't explicitly handle all the possible values that variable can have. For example:

void HandleColinearCase(ColinearPair& curr_minmax, const Point2d& p1, const Point2d& p2, const Point2d& q) {
    switch (getOrder(p1, p2, q)) {
    case Order::Before:
    case Order::EqualSrc:
        UpdateColinearPair(curr_minmax, q, p2);
        break;
    case Order::Between;
        UpdateColinearPair(curr_minmax, p1, p2);
        break;
    case Order::EqualDst:
    case Order::After:
        UpdateColinearPair(curr_minmax, p1, q);
        break;
    }
}

In the above, there is no need for a default case; the compiler can see that you handle all possible results.

Inefficient adding to a convex hull

JDługosz already mentioned that the first argument to ConvexHullAdd() should probably made const. Indeed that looks the case from how it is used:

for (const auto& p : CH1) {
    CH = ConvexHullAdd(CH, p);
}

On the right-hand side of the assignment, p is added to CH if possible, and then a copy of CH is returned. This copy then replaces the original contents of CH. The end result is what you expect, but the copy was made unnecessarily. It would be more efficient to make ConvexHullAdd() void, and to replace the above loop with:

for (const auto& p : CH1) {
    ConvexHullAdd(CH, p);
}

About that tolerance

You hardcoded the tolerance to be 1e-6. This works fine if the spread of your points is much larger than the tolerance. But what if the coordinates are all much smaller than 1e-6? Then isClose() wil always return true, and all expressions of the form a - b > tolerance will be false. The result of getOrder() will then never be Order::Before or Order::Between anymore.

Hardcoding a tolerance value is almost never a good idea. If you have encountered issues (perhaps Cos() is returning infinity or NaN if points are too close), then you should try to solve them in a different way.