m1une's library

This documentation is automatically generated by online-judge-tools/verification-helper

View on GitHub

:heavy_check_mark: Area of Union of Rectangles
(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

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
Back to top page