Area of Union of Rectangles
(geometry/rectangle_union_area.hpp)
- View this file on GitHub
- Last update: 2026-10-05 22:23:07+09:00
- Include:
#include "geometry/rectangle_union_area.hpp"
Overview
rectangle_union_area returns the area covered by at least one axis-aligned
rectangle. Overlaps are counted once. Coordinates may be negative, rectangles
may touch, and zero-width or zero-height rectangles contribute no area.
The implementation sweeps rectangle edges from left to right. A segment tree over compressed y-coordinates maintains the total y-length currently covered. Between consecutive x-events, that length is multiplied by the x-distance.
Interface
template <Coordinate T>
struct AxisAlignedRectangle {
T left;
T right;
T bottom;
T top;
};
template <Coordinate T>
wide_type<T> rectangle_union_area(
const std::vector<AxisAlignedRectangle<T>>& rectangles
);
| Member / Function | Complexity | Description |
|---|---|---|
AxisAlignedRectangle<T>() |
$O(1)$ | Constructs a zero-area rectangle at the origin. |
AxisAlignedRectangle<T>(left, right, bottom, top) |
$O(1)$ | Constructs [left, right] by [bottom, top]. |
rectangle_union_area(rectangles) |
$O(N\log N)$ time and $O(N)$ memory | Returns the union area without modifying the input. |
The constructor coordinates must satisfy left <= right and bottom <= top.
This is checked when the area is calculated. Boundary inclusion does not affect
area; the implementation treats rectangles as half-open during the sweep.
For integral T, wide_type<T> is signed __int128_t, preventing overflow
in coordinate subtraction and multiplication. For floating-point T, the
return type is long double.
Example
#include "geometry/rectangle_union_area.hpp"
#include <iostream>
#include <vector>
int main() {
using Rectangle = m1une::geometry::AxisAlignedRectangle<long long>;
std::vector<Rectangle> rectangles;
rectangles.emplace_back(0, 4, 0, 3);
rectangles.emplace_back(2, 6, 1, 5);
long long area = static_cast<long long>(
m1une::geometry::rectangle_union_area(rectangles)
);
std::cout << area << "\n"; // 24
}
Depends on
Required by
Verified with
verify/geometry/centroid.test.cpp
verify/geometry/geometry_algorithms.test.cpp
verify/geometry/rational.test.cpp
verify/geometry/rectangle_union_area.test.cpp
Code
#ifndef M1UNE_GEOMETRY_RECTANGLE_UNION_AREA_HPP
#define M1UNE_GEOMETRY_RECTANGLE_UNION_AREA_HPP 1
#include "point.hpp"
#include <algorithm>
#include <cassert>
#include <vector>
namespace m1une {
namespace geometry {
template <Coordinate T>
struct AxisAlignedRectangle {
T left;
T right;
T bottom;
T top;
constexpr AxisAlignedRectangle()
: left(0), right(0), bottom(0), top(0) {}
constexpr AxisAlignedRectangle(
T left_value,
T right_value,
T bottom_value,
T top_value
)
: left(left_value),
right(right_value),
bottom(bottom_value),
top(top_value) {}
friend constexpr bool operator==(
const AxisAlignedRectangle&,
const AxisAlignedRectangle&
) = default;
};
namespace rectangle_union_area_detail {
template <Coordinate T>
struct Event {
T x;
T bottom;
T top;
int delta;
Event(T x_value, T bottom_value, T top_value, int delta_value)
: x(x_value),
bottom(bottom_value),
top(top_value),
delta(delta_value) {}
};
template <Coordinate T>
struct CoveredLengthTree {
using Wide = wide_type<T>;
int interval_count;
const std::vector<T>& coordinates;
std::vector<int> cover;
std::vector<Wide> covered;
explicit CoveredLengthTree(const std::vector<T>& values)
: interval_count(std::max(0, int(values.size()) - 1)),
coordinates(values),
cover(std::max(1, 4 * interval_count), 0),
covered(std::max(1, 4 * interval_count), Wide(0)) {}
void pull(int node, int left, int right) {
if (cover[node] > 0) {
covered[node] =
Wide(coordinates[right]) - Wide(coordinates[left]);
} else if (right - left == 1) {
covered[node] = 0;
} else {
covered[node] = covered[2 * node] + covered[2 * node + 1];
}
}
void add(
int query_left,
int query_right,
int delta,
int node,
int left,
int right
) {
if (query_right <= left || right <= query_left) return;
if (query_left <= left && right <= query_right) {
cover[node] += delta;
assert(cover[node] >= 0);
pull(node, left, right);
return;
}
int middle = (left + right) / 2;
add(query_left, query_right, delta, 2 * node, left, middle);
add(
query_left,
query_right,
delta,
2 * node + 1,
middle,
right
);
pull(node, left, right);
}
void add(int left, int right, int delta) {
if (left >= right) return;
assert(0 <= left && left < right && right <= interval_count);
add(left, right, delta, 1, 0, interval_count);
}
Wide length() const {
return interval_count == 0 ? Wide(0) : covered[1];
}
};
} // namespace rectangle_union_area_detail
template <Coordinate T>
wide_type<T> rectangle_union_area(
const std::vector<AxisAlignedRectangle<T>>& rectangles
) {
using Wide = wide_type<T>;
using Event = rectangle_union_area_detail::Event<T>;
std::vector<Event> events;
std::vector<T> y_coordinates;
events.reserve(2 * rectangles.size());
y_coordinates.reserve(2 * rectangles.size());
for (const AxisAlignedRectangle<T>& rectangle : rectangles) {
assert(rectangle.left <= rectangle.right);
assert(rectangle.bottom <= rectangle.top);
if (
rectangle.left == rectangle.right ||
rectangle.bottom == rectangle.top
) {
continue;
}
events.emplace_back(
rectangle.left,
rectangle.bottom,
rectangle.top,
1
);
events.emplace_back(
rectangle.right,
rectangle.bottom,
rectangle.top,
-1
);
y_coordinates.push_back(rectangle.bottom);
y_coordinates.push_back(rectangle.top);
}
if (events.empty()) return Wide(0);
std::sort(y_coordinates.begin(), y_coordinates.end());
y_coordinates.erase(
std::unique(y_coordinates.begin(), y_coordinates.end()),
y_coordinates.end()
);
std::sort(events.begin(), events.end(), [](const Event& a, const Event& b) {
return a.x < b.x;
});
rectangle_union_area_detail::CoveredLengthTree<T> tree(y_coordinates);
Wide area = 0;
T previous_x = events.front().x;
int event_index = 0;
while (event_index < int(events.size())) {
T x = events[event_index].x;
area += tree.length() * (Wide(x) - Wide(previous_x));
int next = event_index;
while (next < int(events.size()) && events[next].x == x) {
int bottom = int(
std::lower_bound(
y_coordinates.begin(),
y_coordinates.end(),
events[next].bottom
) - y_coordinates.begin()
);
int top = int(
std::lower_bound(
y_coordinates.begin(),
y_coordinates.end(),
events[next].top
) - y_coordinates.begin()
);
tree.add(bottom, top, events[next].delta);
next++;
}
previous_x = x;
event_index = next;
}
return area;
}
} // namespace geometry
} // namespace m1une
#endif // M1UNE_GEOMETRY_RECTANGLE_UNION_AREA_HPP#line 1 "geometry/rectangle_union_area.hpp"
#line 1 "geometry/point.hpp"
#include <cmath>
#include <concepts>
#include <cassert>
#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 5 "geometry/rectangle_union_area.hpp"
#include <algorithm>
#line 8 "geometry/rectangle_union_area.hpp"
#include <vector>
namespace m1une {
namespace geometry {
template <Coordinate T>
struct AxisAlignedRectangle {
T left;
T right;
T bottom;
T top;
constexpr AxisAlignedRectangle()
: left(0), right(0), bottom(0), top(0) {}
constexpr AxisAlignedRectangle(
T left_value,
T right_value,
T bottom_value,
T top_value
)
: left(left_value),
right(right_value),
bottom(bottom_value),
top(top_value) {}
friend constexpr bool operator==(
const AxisAlignedRectangle&,
const AxisAlignedRectangle&
) = default;
};
namespace rectangle_union_area_detail {
template <Coordinate T>
struct Event {
T x;
T bottom;
T top;
int delta;
Event(T x_value, T bottom_value, T top_value, int delta_value)
: x(x_value),
bottom(bottom_value),
top(top_value),
delta(delta_value) {}
};
template <Coordinate T>
struct CoveredLengthTree {
using Wide = wide_type<T>;
int interval_count;
const std::vector<T>& coordinates;
std::vector<int> cover;
std::vector<Wide> covered;
explicit CoveredLengthTree(const std::vector<T>& values)
: interval_count(std::max(0, int(values.size()) - 1)),
coordinates(values),
cover(std::max(1, 4 * interval_count), 0),
covered(std::max(1, 4 * interval_count), Wide(0)) {}
void pull(int node, int left, int right) {
if (cover[node] > 0) {
covered[node] =
Wide(coordinates[right]) - Wide(coordinates[left]);
} else if (right - left == 1) {
covered[node] = 0;
} else {
covered[node] = covered[2 * node] + covered[2 * node + 1];
}
}
void add(
int query_left,
int query_right,
int delta,
int node,
int left,
int right
) {
if (query_right <= left || right <= query_left) return;
if (query_left <= left && right <= query_right) {
cover[node] += delta;
assert(cover[node] >= 0);
pull(node, left, right);
return;
}
int middle = (left + right) / 2;
add(query_left, query_right, delta, 2 * node, left, middle);
add(
query_left,
query_right,
delta,
2 * node + 1,
middle,
right
);
pull(node, left, right);
}
void add(int left, int right, int delta) {
if (left >= right) return;
assert(0 <= left && left < right && right <= interval_count);
add(left, right, delta, 1, 0, interval_count);
}
Wide length() const {
return interval_count == 0 ? Wide(0) : covered[1];
}
};
} // namespace rectangle_union_area_detail
template <Coordinate T>
wide_type<T> rectangle_union_area(
const std::vector<AxisAlignedRectangle<T>>& rectangles
) {
using Wide = wide_type<T>;
using Event = rectangle_union_area_detail::Event<T>;
std::vector<Event> events;
std::vector<T> y_coordinates;
events.reserve(2 * rectangles.size());
y_coordinates.reserve(2 * rectangles.size());
for (const AxisAlignedRectangle<T>& rectangle : rectangles) {
assert(rectangle.left <= rectangle.right);
assert(rectangle.bottom <= rectangle.top);
if (
rectangle.left == rectangle.right ||
rectangle.bottom == rectangle.top
) {
continue;
}
events.emplace_back(
rectangle.left,
rectangle.bottom,
rectangle.top,
1
);
events.emplace_back(
rectangle.right,
rectangle.bottom,
rectangle.top,
-1
);
y_coordinates.push_back(rectangle.bottom);
y_coordinates.push_back(rectangle.top);
}
if (events.empty()) return Wide(0);
std::sort(y_coordinates.begin(), y_coordinates.end());
y_coordinates.erase(
std::unique(y_coordinates.begin(), y_coordinates.end()),
y_coordinates.end()
);
std::sort(events.begin(), events.end(), [](const Event& a, const Event& b) {
return a.x < b.x;
});
rectangle_union_area_detail::CoveredLengthTree<T> tree(y_coordinates);
Wide area = 0;
T previous_x = events.front().x;
int event_index = 0;
while (event_index < int(events.size())) {
T x = events[event_index].x;
area += tree.length() * (Wide(x) - Wide(previous_x));
int next = event_index;
while (next < int(events.size()) && events[next].x == x) {
int bottom = int(
std::lower_bound(
y_coordinates.begin(),
y_coordinates.end(),
events[next].bottom
) - y_coordinates.begin()
);
int top = int(
std::lower_bound(
y_coordinates.begin(),
y_coordinates.end(),
events[next].top
) - y_coordinates.begin()
);
tree.add(bottom, top, events[next].delta);
next++;
}
previous_x = x;
event_index = next;
}
return area;
}
} // namespace geometry
} // namespace m1une