Minkowski Sum
(geometry/minkowski_sum.hpp)
- View this file on GitHub
- Last update: 2026-10-05 22:23:07+09:00
- Include:
#include "geometry/minkowski_sum.hpp"
Overview
minkowski_sum constructs the Minkowski sum of two ordered convex polygons:
the set of all points a + b for a point a in the first polygon and a point
b in the second polygon.
The input boundaries may be clockwise or counterclockwise and may contain a repeated closing point, consecutive duplicates, or redundant collinear vertices. Points and segments are supported as degenerate convex polygons.
Function
template <Coordinate T>
std::vector<Point<T>> minkowski_sum(
std::vector<Point<T>> first,
std::vector<Point<T>> second,
long double eps = 1e-12L
);
| Function | Description | Complexity |
|---|---|---|
minkowski_sum(first, second, eps) |
Returns the normalized boundary of the Minkowski sum. | $O(N+M)$ time and memory |
Both inputs must be nonempty and must describe convex boundaries in cyclic
order. The result is counterclockwise, starts at its lowest (y, x) vertex,
has no repeated closing point, and omits redundant collinear vertices.
The return type keeps the input coordinate type. Coordinate addition and edge
subtraction must therefore fit T. Geometric predicates use wide_type<T>;
integral coordinates consequently use signed 128-bit cross products. For
floating-point coordinates, eps controls normalization of collinear points.
Example
#include "geometry/minkowski_sum.hpp"
#include <iostream>
#include <vector>
int main() {
using Point = m1une::geometry::Point<long long>;
std::vector<Point> square;
square.emplace_back(0, 0);
square.emplace_back(2, 0);
square.emplace_back(2, 2);
square.emplace_back(0, 2);
std::vector<Point> segment;
segment.emplace_back(0, 0);
segment.emplace_back(3, 0);
const auto sum = m1une::geometry::minkowski_sum(square, segment);
for (const Point& point : sum) {
std::cout << point.x << " " << point.y << "\n";
}
}
Depends on
- geometry/detail/convex_polygon_normalize.hpp
- geometry/detail/floating_predicate.hpp
- 2D Point and Predicates
(geometry/point.hpp)
Required by
Verified with
- verify/geometry/centroid.test.cpp
- verify/geometry/convex_decomposition.test.cpp
- verify/geometry/convex_diameter.test.cpp
- verify/geometry/convex_polygon.test.cpp
- verify/geometry/geometry_algorithms.test.cpp
- verify/geometry/is_convex_polygon.test.cpp
- verify/geometry/minkowski_sum.test.cpp
- verify/geometry/polygon_operations.test.cpp
- verify/geometry/rational.test.cpp
- verify/geometry/steiner_convex_decomposition.test.cpp
Code
#ifndef M1UNE_GEOMETRY_MINKOWSKI_SUM_HPP
#define M1UNE_GEOMETRY_MINKOWSKI_SUM_HPP 1
#include <cassert>
#include <cstddef>
#include <utility>
#include <vector>
#include "detail/convex_polygon_normalize.hpp"
namespace m1une {
namespace geometry {
// Returns the normalized boundary of the Minkowski sum of two nonempty
// ordered convex polygons.
template <Coordinate T>
std::vector<Point<T>> minkowski_sum(
std::vector<Point<T>> first,
std::vector<Point<T>> second,
long double eps = 1e-12L
) {
assert(!first.empty());
assert(!second.empty());
first = convex_polygon_detail::normalize_convex_boundary(
std::move(first),
eps
);
second = convex_polygon_detail::normalize_convex_boundary(
std::move(second),
eps
);
if (first.size() == 1 || second.size() == 1) {
if (second.size() == 1) std::swap(first, second);
for (Point<T>& point : second) point += first[0];
return convex_polygon_detail::normalize_convex_boundary(
std::move(second),
eps
);
}
std::vector<Point<T>> first_edges;
std::vector<Point<T>> second_edges;
first_edges.reserve(first.size());
second_edges.reserve(second.size());
for (std::size_t index = 0; index < first.size(); ++index) {
first_edges.push_back(
first[(index + 1) % first.size()] - first[index]
);
}
for (std::size_t index = 0; index < second.size(); ++index) {
second_edges.push_back(
second[(index + 1) % second.size()] - second[index]
);
}
Point<T> current = first.front() + second.front();
std::vector<Point<T>> result;
result.reserve(first.size() + second.size());
result.push_back(current);
std::size_t first_index = 0;
std::size_t second_index = 0;
while (
first_index < first_edges.size() ||
second_index < second_edges.size()
) {
Point<T> step;
if (first_index == first_edges.size()) {
step = second_edges[second_index++];
} else if (second_index == second_edges.size()) {
step = first_edges[first_index++];
} else {
const auto turn = cross(
first_edges[first_index],
second_edges[second_index]
);
if (turn > 0) {
step = first_edges[first_index++];
} else if (turn < 0) {
step = second_edges[second_index++];
} else {
step = first_edges[first_index++] +
second_edges[second_index++];
}
}
current += step;
if (
first_index < first_edges.size() ||
second_index < second_edges.size()
) {
result.push_back(current);
}
}
return convex_polygon_detail::normalize_convex_boundary(
std::move(result),
eps
);
}
} // namespace geometry
} // namespace m1une
#endif // M1UNE_GEOMETRY_MINKOWSKI_SUM_HPP#line 1 "geometry/minkowski_sum.hpp"
#include <cassert>
#include <cstddef>
#include <utility>
#include <vector>
#line 1 "geometry/detail/convex_polygon_normalize.hpp"
#include <algorithm>
#line 8 "geometry/detail/convex_polygon_normalize.hpp"
#line 1 "geometry/point.hpp"
#include <cmath>
#include <concepts>
#line 7 "geometry/point.hpp"
#include <type_traits>
#line 1 "geometry/detail/floating_predicate.hpp"
namespace m1une {
namespace geometry {
namespace predicate_detail {
template <typename T>
constexpr T absolute(T value) {
return value < T(0) ? -value : value;
}
template <typename T>
constexpr T max_value(T first, T second) {
return first < second ? second : first;
}
template <typename T>
constexpr T vector_scale(T x, T y) {
return max_value(absolute(x), absolute(y));
}
template <bool Exact, typename T>
constexpr int scaled_sign(T value, T scale, long double eps) {
if constexpr (Exact) {
return (value > T(0)) - (value < T(0));
} else {
const T tolerance = T(eps) * scale;
return (value > tolerance) - (value < -tolerance);
}
}
template <bool Exact, typename T>
constexpr T determinant_scale(T ax, T ay, T bx, T by) {
if constexpr (Exact) {
return T(0);
} else {
return vector_scale(ax, ay) * vector_scale(bx, by);
}
}
template <bool Exact, typename T>
constexpr int determinant_sign(
T ax,
T ay,
T bx,
T by,
long double eps
) {
const T determinant = ax * by - ay * bx;
return scaled_sign<Exact>(
determinant,
determinant_scale<Exact>(ax, ay, bx, by),
eps
);
}
template <bool Exact, typename T>
constexpr int orientation_sign(
T direction_x,
T direction_y,
T offset_x,
T offset_y,
long double eps
) {
const T determinant =
direction_x * offset_y - direction_y * offset_x;
T scale = T(0);
if constexpr (!Exact) {
const T direction_scale =
vector_scale(direction_x, direction_y);
scale = direction_scale * max_value(
direction_scale,
vector_scale(offset_x, offset_y)
);
}
return scaled_sign<Exact>(determinant, scale, eps);
}
template <bool Exact, typename T>
constexpr int dot_sign(
T ax,
T ay,
T bx,
T by,
long double eps
) {
const T value = ax * bx + ay * by;
T scale = T(0);
if constexpr (!Exact) {
scale = vector_scale(ax, ay) * vector_scale(bx, by);
}
return scaled_sign<Exact>(value, scale, eps);
}
} // namespace predicate_detail
} // namespace geometry
} // namespace m1une
#line 10 "geometry/point.hpp"
namespace m1une {
namespace geometry {
template <typename T>
concept Coordinate = !std::same_as<std::remove_cv_t<T>, bool> &&
(std::is_arithmetic_v<T> ||
(std::copyable<T> && std::totally_ordered<T> && requires(T a, T b) {
T(0);
T(1);
static_cast<long double>(a);
{ +a } -> std::same_as<T>;
{ -a } -> std::same_as<T>;
{ a + b } -> std::same_as<T>;
{ a - b } -> std::same_as<T>;
{ a * b } -> std::same_as<T>;
{ a / b } -> std::same_as<T>;
{ a += b } -> std::same_as<T&>;
{ a -= b } -> std::same_as<T&>;
}));
// Custom coordinate types keep their own exact arithmetic.
template <typename T>
concept ExactCoordinate = Coordinate<T> && !std::floating_point<T>;
template <Coordinate T>
using wide_type = std::conditional_t<std::integral<T>, __int128_t,
std::conditional_t<std::floating_point<T>, long double, T>>;
template <Coordinate T>
struct Point {
T x;
T y;
constexpr Point() : x(0), y(0) {}
constexpr Point(T x_value, T y_value) : x(x_value), y(y_value) {}
template <Coordinate U>
explicit constexpr Point(const Point<U>& other)
: x(static_cast<T>(other.x)), y(static_cast<T>(other.y)) {}
constexpr Point& operator+=(const Point& other) {
x += other.x;
y += other.y;
return *this;
}
constexpr Point& operator-=(const Point& other) {
x -= other.x;
y -= other.y;
return *this;
}
constexpr Point operator+() const {
return *this;
}
constexpr Point operator-() const {
return Point(-x, -y);
}
friend constexpr Point operator+(Point left, const Point& right) {
return left += right;
}
friend constexpr Point operator-(Point left, const Point& right) {
return left -= right;
}
friend constexpr bool operator==(const Point&, const Point&) = default;
friend constexpr bool operator<(const Point& left, const Point& right) {
if (left.x != right.x) return left.x < right.x;
return left.y < right.y;
}
};
template <Coordinate T>
constexpr Point<long double> centroid(const Point<T>& point) {
return Point<long double>(point);
}
template <Coordinate T, typename Scalar>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator*(const Point<T>& point, Scalar scalar) {
using Result = std::common_type_t<T, Scalar>;
return Point<Result>(
Result(point.x) * Result(scalar),
Result(point.y) * Result(scalar)
);
}
template <typename Scalar, Coordinate T>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator*(Scalar scalar, const Point<T>& point) {
return point * scalar;
}
template <Coordinate T, typename Scalar>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator/(const Point<T>& point, Scalar scalar) {
using Result = std::common_type_t<T, Scalar>;
return Point<Result>(
Result(point.x) / Result(scalar),
Result(point.y) / Result(scalar)
);
}
template <Coordinate T>
constexpr wide_type<T> dot(const Point<T>& a, const Point<T>& b) {
using W = wide_type<T>;
return W(a.x) * W(b.x) + W(a.y) * W(b.y);
}
template <Coordinate T>
constexpr wide_type<T> cross(const Point<T>& a, const Point<T>& b) {
using W = wide_type<T>;
return W(a.x) * W(b.y) - W(a.y) * W(b.x);
}
template <Coordinate T>
constexpr wide_type<T> cross(
const Point<T>& origin,
const Point<T>& a,
const Point<T>& b
) {
using W = wide_type<T>;
W ax = W(a.x) - W(origin.x);
W ay = W(a.y) - W(origin.y);
W bx = W(b.x) - W(origin.x);
W by = W(b.y) - W(origin.y);
return ax * by - ay * bx;
}
template <Coordinate T>
constexpr wide_type<T> norm2(const Point<T>& point) {
return dot(point, point);
}
template <Coordinate T>
constexpr wide_type<T> distance2(const Point<T>& a, const Point<T>& b) {
using W = wide_type<T>;
W dx = W(a.x) - W(b.x);
W dy = W(a.y) - W(b.y);
return dx * dx + dy * dy;
}
template <Coordinate T>
long double norm(const Point<T>& point) {
return std::hypot(
static_cast<long double>(point.x),
static_cast<long double>(point.y)
);
}
template <Coordinate T>
long double distance(const Point<T>& a, const Point<T>& b) {
return std::hypot(
static_cast<long double>(a.x) - static_cast<long double>(b.x),
static_cast<long double>(a.y) - static_cast<long double>(b.y)
);
}
template <Coordinate T, typename M, typename N>
requires (std::is_arithmetic_v<M> || Coordinate<M>) &&
(std::is_arithmetic_v<N> || Coordinate<N>)
constexpr Point<long double> internal_division_point(
const Point<T>& a,
const Point<T>& b,
M m,
N n
) {
long double first_ratio = static_cast<long double>(m);
long double second_ratio = static_cast<long double>(n);
long double denominator = first_ratio + second_ratio;
assert(denominator != 0);
Point<long double> first(a);
Point<long double> direction = Point<long double>(b) - first;
return first + direction * (first_ratio / denominator);
}
template <Coordinate T, typename M, typename N>
requires (std::is_arithmetic_v<M> || Coordinate<M>) &&
(std::is_arithmetic_v<N> || Coordinate<N>)
constexpr Point<long double> external_division_point(
const Point<T>& a,
const Point<T>& b,
M m,
N n
) {
long double first_ratio = static_cast<long double>(m);
long double second_ratio = static_cast<long double>(n);
long double denominator = first_ratio - second_ratio;
assert(denominator != 0);
Point<long double> first(a);
Point<long double> direction = Point<long double>(b) - first;
return first + direction * (first_ratio / denominator);
}
template <Coordinate T>
constexpr int sign(wide_type<T> value, long double eps = 1e-12L) {
return predicate_detail::scaled_sign<ExactCoordinate<T>>(
value,
wide_type<T>(1),
eps
);
}
template <Coordinate T>
constexpr int orientation(
const Point<T>& a,
const Point<T>& b,
const Point<T>& c,
long double eps = 1e-12L
) {
using W = wide_type<T>;
const W first_x = W(b.x) - W(a.x);
const W first_y = W(b.y) - W(a.y);
const W second_x = W(c.x) - W(a.x);
const W second_y = W(c.y) - W(a.y);
return predicate_detail::orientation_sign<ExactCoordinate<T>>(
first_x,
first_y,
second_x,
second_y,
eps
);
}
template <Coordinate T>
constexpr bool collinear(
const Point<T>& a,
const Point<T>& b,
const Point<T>& c,
long double eps = 1e-12L
) {
return orientation(a, b, c, eps) == 0;
}
template <Coordinate T>
Point<long double> rotate(const Point<T>& point, long double angle) {
long double cosine = std::cos(angle);
long double sine = std::sin(angle);
return Point<long double>(
static_cast<long double>(point.x) * cosine -
static_cast<long double>(point.y) * sine,
static_cast<long double>(point.x) * sine +
static_cast<long double>(point.y) * cosine
);
}
template <Coordinate T>
Point<long double> normalized(const Point<T>& point) {
long double length = norm(point);
assert(length != 0);
return Point<long double>(
static_cast<long double>(point.x) / length,
static_cast<long double>(point.y) / length
);
}
} // namespace geometry
} // namespace m1une
#line 10 "geometry/detail/convex_polygon_normalize.hpp"
namespace m1une {
namespace geometry {
namespace convex_polygon_detail {
template <Coordinate T>
wide_type<T> boundary_area2(const std::vector<Point<T>>& polygon) {
wide_type<T> result = 0;
for (std::size_t index = 0; index < polygon.size(); ++index) {
result += cross(
polygon[index],
polygon[(index + 1) % polygon.size()]
);
}
return result;
}
template <Coordinate T>
std::vector<Point<T>> normalize_convex_boundary(
std::vector<Point<T>> polygon,
long double eps
) {
if (polygon.size() >= 2 && polygon.front() == polygon.back()) {
polygon.pop_back();
}
polygon.erase(
std::unique(polygon.begin(), polygon.end()),
polygon.end()
);
if (polygon.size() >= 2 && polygon.front() == polygon.back()) {
polygon.pop_back();
}
if (polygon.size() <= 1) return polygon;
if (
polygon.size() >= 3 &&
sign<T>(boundary_area2(polygon), eps) < 0
) {
std::reverse(polygon.begin(), polygon.end());
}
const auto start = std::min_element(
polygon.begin(),
polygon.end(),
[](const Point<T>& first, const Point<T>& second) {
if (first.y != second.y) return first.y < second.y;
return first.x < second.x;
}
);
std::rotate(polygon.begin(), start, polygon.end());
if (polygon.size() >= 3) {
std::vector<Point<T>> cleaned;
const std::size_t size = polygon.size();
cleaned.reserve(size);
for (std::size_t index = 0; index < size; ++index) {
const Point<T>& previous = polygon[(index + size - 1) % size];
const Point<T>& current = polygon[index];
const Point<T>& next = polygon[(index + 1) % size];
if (
orientation(previous, current, next, eps) != 0 ||
sign<T>(dot(current - previous, next - current), eps) < 0
) {
cleaned.push_back(current);
}
}
polygon = std::move(cleaned);
}
return polygon;
}
} // namespace convex_polygon_detail
} // namespace geometry
} // namespace m1une
#line 10 "geometry/minkowski_sum.hpp"
namespace m1une {
namespace geometry {
// Returns the normalized boundary of the Minkowski sum of two nonempty
// ordered convex polygons.
template <Coordinate T>
std::vector<Point<T>> minkowski_sum(
std::vector<Point<T>> first,
std::vector<Point<T>> second,
long double eps = 1e-12L
) {
assert(!first.empty());
assert(!second.empty());
first = convex_polygon_detail::normalize_convex_boundary(
std::move(first),
eps
);
second = convex_polygon_detail::normalize_convex_boundary(
std::move(second),
eps
);
if (first.size() == 1 || second.size() == 1) {
if (second.size() == 1) std::swap(first, second);
for (Point<T>& point : second) point += first[0];
return convex_polygon_detail::normalize_convex_boundary(
std::move(second),
eps
);
}
std::vector<Point<T>> first_edges;
std::vector<Point<T>> second_edges;
first_edges.reserve(first.size());
second_edges.reserve(second.size());
for (std::size_t index = 0; index < first.size(); ++index) {
first_edges.push_back(
first[(index + 1) % first.size()] - first[index]
);
}
for (std::size_t index = 0; index < second.size(); ++index) {
second_edges.push_back(
second[(index + 1) % second.size()] - second[index]
);
}
Point<T> current = first.front() + second.front();
std::vector<Point<T>> result;
result.reserve(first.size() + second.size());
result.push_back(current);
std::size_t first_index = 0;
std::size_t second_index = 0;
while (
first_index < first_edges.size() ||
second_index < second_edges.size()
) {
Point<T> step;
if (first_index == first_edges.size()) {
step = second_edges[second_index++];
} else if (second_index == second_edges.size()) {
step = first_edges[first_index++];
} else {
const auto turn = cross(
first_edges[first_index],
second_edges[second_index]
);
if (turn > 0) {
step = first_edges[first_index++];
} else if (turn < 0) {
step = second_edges[second_index++];
} else {
step = first_edges[first_index++] +
second_edges[second_index++];
}
}
current += step;
if (
first_index < first_edges.size() ||
second_index < second_edges.size()
) {
result.push_back(current);
}
}
return convex_polygon_detail::normalize_convex_boundary(
std::move(result),
eps
);
}
} // namespace geometry
} // namespace m1une