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';
 }
}
````

• JDługosz: Are you required to use a specific (non-current) version of C++? That is, if you can't use C++20 features we won't point out that you're doing things the old (hard, verbose) way.
• frozenca: @JDługosz No, the opposite is true: I prefer to use the most recent version of C++, because the codes I'm posting here is purely for fun and hobby projects, not for production code that should consider backward compatibility. If you see outdated practices in my code, then the reason is simply I don't know the best up-to-date practice. Pointing out that would be welcomed.
• JDługosz: see en.cppreference.com/w/cpp/language/default_comparisons for a welcome improvement! See Number 4 on en.cppreference.com/w/cpp/language/…

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.

• Toby Speight: Your first point is interesting, because I go the opposite way - include the project-specific libraries first, to help spot any hidden dependencies they may have.
• indi: @TobySpeight That makes sense in testing code (and that’s how I do it: header-under-test first, usually), but it causes problems generally. Including std headers first, and in a specific order (usually alphabetical, with some tweaks sometimes (like maybe <version> first)) makes for better build consistency, and makes it easier to cache/precompile std headers.
• indi: "Well, your constructors are not constexpr-friendly. There's no reason why it can't be!" Rather than making the constructors constexpr, it would be better to just delete them. They’re unnecessary.
• JDługosz: @indi hmm, is an implicitly declared (or defaulted on first declaration) default constructor going to be constexpr? Yes, it is. But no such wording is present for the aggregate initializer, so the two-argument constructor (and by extension a vector of them that's all literal) won't be done at compile time.
• JDługosz: hmm again... I see en.cppreference.com/w/cpp/named_req/LiteralType indicates that it may be an aggregate. So the aggregate initializer would be allowed in a variable declared as constexpr.

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.