Given is a set S of n axis-aligned cubes. The task is to find the volume of the union of all cubes in S. This means that every volume-overlap of 2 or more cubes is only counted once. The set specifically contains all the coordinates of all the cubes.
I have found several papers regarding the subject, presenting algorithms to complete the task. This paper for example generalizes the problem to any dimension d rather than the trivial d=3, and to boxes rather than cubes. This paper as well as a few others solve the problem in time O(n^1.5) or slightly better. Another paper which looks promising and is specific to 3d-cubes is this one which solves the task in O(n^4/3 log n). But the papers seem rather complex, at least for me, and I cannot follow them clearly.
Is there any relatively simple pseudocode or procedure that I can follow to implement this idea? I am looking for a set of instructions, what exactly to do with the cubes. Any implementation will be also excellent. And O(n^2) or even O(n^3) are fine.
Currently, my approach was to compute the total volume, i.e the sum of all volumes of all the cubes, and then compute the overlap of every 2 cubes, and reduce it from the total volume. The problem is that every such overlap may (or may not) be of a different pair of cubes, meaning an overlap can be for example shared by 5 cubes. In this approach the overlap will be counted 10 times rather just once. So I was thinking maybe an inclusion-exclusion principle may prove itself useful, but I don't know exactly how it may be implemented specifically. Computing every overlap (naively) is already O(n^2), but doable. Is there any better way to compute all such overlaps?
Anyways, I don't assume this is a useful approach, to continue along these lines.
I implemented Bentley's algorithm (O(n^2 log n)) in C++. (I know you wanted Java, but C++ is my main ax at work, and templates are just too useful here given that I was considering working my way up to Overmars and Yap.)
// Import some basic stuff from the standard library.
#include <algorithm>
#include <cassert>
#include <cmath>
#include <iostream>
#include <limits>
#include <memory>
#include <optional>
#include <random>
#include <tuple>
#include <utility>
#include <vector>
// Define vocabulary types.
class Interval {
public:
Interval(double a, double b) : min_(std::fmin(a, b)), max_(std::fmax(a, b)) {}
double min() const { return min_; }
double max() const { return max_; }
private:
double min_, max_;
};
// Cartesian product of an interval and a set.
template <typename Projection>
class IntervalTimes {
public:
IntervalTimes(Interval interval, Projection projection)
: interval_(interval), projection_(projection) {}
Interval interval() const { return interval_; }
Projection projection() const { return projection_; }
private:
Interval interval_;
Projection projection_;
};
// Zero-dimensional base case.
struct Interval0 {};
using Interval1 = IntervalTimes<Interval0>;
using Interval2 = IntervalTimes<Interval1>;
using Interval3 = IntervalTimes<Interval2>;
using Box = Interval3;
Box MakeBox(Interval x, Interval y, Interval z) {
return IntervalTimes{x, IntervalTimes{y, IntervalTimes{z, Interval0{}}}};
}
// Define basic operations.
double Length(Interval interval) { return interval.max() - interval.min(); }
double Measure(Interval0) { return 1; }
template <typename Projection>
double Measure(IntervalTimes<Projection> product) {
return Length(product.interval()) * Measure(product.projection());
}
bool Contains(Interval interval, double x) {
return interval.min() < x && x < interval.max();
}
bool Contains(Interval i1, Interval i2) {
return i1.min() <= i2.min() && i2.max() <= i1.max();
}
bool Contains(Box box, double x, double y, double z) {
return Contains(box.interval(), x) &&
Contains(box.projection().interval(), y) &&
Contains(box.projection().projection().interval(), z);
}
bool Intersects(Interval i1, Interval i2) {
return std::fmax(i1.min(), i2.min()) < std::fmin(i1.max(), i2.max());
}
template <typename Projection>
std::vector<Projection> ProjectEach(
const std::vector<IntervalTimes<Projection>>& objects) {
std::vector<Projection> projections;
projections.reserve(objects.size());
for (const IntervalTimes<Projection>& object : objects) {
projections.push_back(object.projection());
}
return projections;
}
template <typename T>
std::vector<T> Select(const std::vector<bool>& included,
const std::vector<T>& objects) {
std::vector<T> selected;
for (std::size_t j = 0; j < included.size() && j < objects.size(); j++) {
if (included[j]) selected.push_back(objects[j]);
}
return selected;
}
// Returns the unique x values that appear in objects, sorted.
template <typename Projection>
std::vector<double> Boundaries(
const std::vector<IntervalTimes<Projection>>& objects) {
std::vector<double> boundaries;
boundaries.reserve(objects.size() * 2);
for (const IntervalTimes<Projection>& object : objects) {
boundaries.push_back(object.interval().min());
boundaries.push_back(object.interval().max());
}
std::sort(boundaries.begin(), boundaries.end());
boundaries.erase(std::unique(boundaries.begin(), boundaries.end()),
boundaries.end());
return boundaries;
}
// The basic offline algorithm for d dimensions uses the online algorithm for
// d-1 dimensions. Each object gives rise to two events. We sweep over the
// events, integrating as we go using the online algorithm.
template <typename Object>
class OnlineMeasure {
public:
virtual ~OnlineMeasure() = default;
virtual void Initialize(std::vector<Object>) {}
// Adds the object at index j in the objects presented to Initialize().
virtual void Add(std::size_t j) = 0;
// Removes the object at index j in the objects presented to Initialize().
virtual void Remove(std::size_t j) = 0;
// Returns the measure of the union of the objects added but not removed.
virtual double Measure() const = 0;
};
enum class Side { kMin, kMax };
// {x, side, index}.
using Event = std::tuple<double, Side, std::size_t>;
template <typename Projection>
double OfflineMeasure(const std::vector<IntervalTimes<Projection>>& objects,
OnlineMeasure<Projection>& online_measure) {
// Construct the events and sort them by x with min before max.
std::vector<Event> events;
events.reserve(objects.size() * 2);
for (std::size_t j = 0; j < objects.size(); j++) {
Interval interval = objects[j].interval();
events.push_back({interval.min(), Side::kMin, j});
events.push_back({interval.max(), Side::kMax, j});
}
std::sort(events.begin(), events.end());
// Sweep x to integrate.
double measure = 0;
std::optional<double> previous_x;
online_measure.Initialize(ProjectEach(objects));
for (const auto& [x, side, j] : events) {
if (previous_x) measure += (x - *previous_x) * online_measure.Measure();
previous_x = x;
switch (side) {
case Side::kMin:
online_measure.Add(j);
break;
case Side::kMax:
online_measure.Remove(j);
break;
}
}
return measure;
}
// We can use the algorithm above as a slow online algorithm.
template <typename Projection>
class OfflineOnlineMeasure : public OnlineMeasure<IntervalTimes<Projection>> {
public:
OfflineOnlineMeasure(
std::unique_ptr<OnlineMeasure<Projection>> online_measure)
: online_measure_(std::move(online_measure)) {}
void Initialize(std::vector<IntervalTimes<Projection>> objects) override {
objects_ = std::move(objects);
included_.assign(objects_.size(), false);
}
void Add(std::size_t j) override { included_.at(j) = true; }
void Remove(std::size_t j) override { included_.at(j) = false; }
double Measure() const override {
return OfflineMeasure(Select(included_, objects_), *online_measure_);
}
private:
std::unique_ptr<OnlineMeasure<Projection>> online_measure_;
std::vector<bool> included_;
std::vector<IntervalTimes<Projection>> objects_;
};
// Zero-dimensional base case for Klee's algorithm.
class KleeOnlineMeasure : public OnlineMeasure<Interval0> {
public:
void Add(std::size_t) override { multiplicity_++; }
void Remove(std::size_t) override { multiplicity_--; }
double Measure() const override { return multiplicity_ > 0 ? 1 : 0; }
private:
std::size_t multiplicity_ = 0;
};
double KleeMeasure(const std::vector<Box>& boxes) {
std::unique_ptr<OnlineMeasure<Interval0>> measure0 =
std::make_unique<KleeOnlineMeasure>();
std::unique_ptr<OnlineMeasure<Interval1>> measure1 =
std::make_unique<OfflineOnlineMeasure<Interval0>>(std::move(measure0));
OfflineOnlineMeasure<Interval1> measure2(std::move(measure1));
return OfflineMeasure(boxes, measure2);
}
// The fundamental insight into Bentley's algorithm is a segment tree that
// solves the online problem in one dimension.
class Segment {
public:
explicit Segment(Interval interval)
: left_{nullptr},
right_{nullptr},
interval_{interval},
multiplicity_{0},
descendant_length_{0} {}
Segment(std::unique_ptr<Segment> left, std::unique_ptr<Segment> right)
: left_{std::move(left)},
right_{std::move(right)},
interval_{left_->interval_.min(), right_->interval_.max()},
multiplicity_{0},
descendant_length_{left_->LengthOfUnion() + right_->LengthOfUnion()} {
assert(left_->interval_.max() == right_->interval_.min());
}
double LengthOfUnion() const {
return multiplicity_ > 0 ? Length(interval_) : descendant_length_;
}
void Add(const Interval i) { Increment(i, 1); }
void Remove(const Interval i) { Increment(i, -1); }
private:
void Increment(const Interval i, std::size_t delta) {
if (Contains(i, interval_)) {
multiplicity_ += delta;
} else if (Intersects(i, interval_)) {
left_->Increment(i, delta);
right_->Increment(i, delta);
descendant_length_ = left_->LengthOfUnion() + right_->LengthOfUnion();
}
}
// Children.
std::unique_ptr<Segment> left_;
std::unique_ptr<Segment> right_;
// Interval.
Interval interval_;
// Adds minus removes for this whole segment.
std::size_t multiplicity_;
// Total length from proper descendants.
double descendant_length_;
};
std::unique_ptr<Segment> ConstructSegmentTree(
const std::vector<double>& boundaries) {
if (boundaries.empty()) return nullptr;
std::vector<std::unique_ptr<Segment>> segments;
segments.reserve(boundaries.size() - 1);
for (std::size_t j = 1; j < boundaries.size(); j++) {
segments.push_back(
std::make_unique<Segment>(Interval{boundaries[j - 1], boundaries[j]}));
}
while (segments.size() > 1) {
std::vector<std::unique_ptr<Segment>> parent_segments;
parent_segments.reserve(segments.size() / 2 + segments.size() % 2);
for (std::size_t j = 1; j < segments.size(); j += 2) {
parent_segments.push_back(std::make_unique<Segment>(
std::move(segments[j - 1]), std::move(segments[j])));
}
if (segments.size() % 2 == 1) {
parent_segments.push_back(std::move(segments.back()));
}
segments = std::move(parent_segments);
}
return std::move(segments.front());
}
class BentleyOnlineMeasure : public OnlineMeasure<Interval1> {
public:
void Initialize(std::vector<Interval1> intervals) override {
intervals_ = std::move(intervals);
root_ = ConstructSegmentTree(Boundaries(intervals_));
}
void Add(std::size_t j) override { root_->Add(intervals_.at(j).interval()); }
void Remove(std::size_t j) override {
root_->Remove(intervals_.at(j).interval());
}
double Measure() const override {
return root_ != nullptr ? root_->LengthOfUnion() : 0;
}
private:
std::vector<Interval1> intervals_;
std::unique_ptr<Segment> root_;
};
double BentleyMeasure(const std::vector<Box>& boxes) {
std::unique_ptr<OnlineMeasure<Interval1>> measure1 =
std::make_unique<BentleyOnlineMeasure>();
OfflineOnlineMeasure<Interval1> measure2(std::move(measure1));
return OfflineMeasure(boxes, measure2);
}
int main() {
std::default_random_engine gen;
std::uniform_real_distribution<double> uniform(0, 1);
std::vector<Box> boxes;
static constexpr std::size_t kBoxCount = 20;
boxes.reserve(kBoxCount);
for (std::size_t j = 0; j < kBoxCount; j++) {
boxes.push_back(MakeBox({uniform(gen), uniform(gen)},
{uniform(gen), uniform(gen)},
{uniform(gen), uniform(gen)}));
}
std::cout << KleeMeasure(boxes) << "\n";
std::cout << BentleyMeasure(boxes) << "\n";
double hits = 0;
static constexpr std::size_t kSampleCount = 1000000;
for (std::size_t j = 0; j < kSampleCount; j++) {
const double x = uniform(gen);
const double y = uniform(gen);
const double z = uniform(gen);
for (const Box& box : boxes) {
if (Contains(box, x, y, z)) {
hits++;
break;
}
}
}
std::cout << hits / kSampleCount << "\n";
}