计算轴对齐立方体并集的体积

Question

给定一组轴对齐的立方体S。任务是求出所有立方体的并n集体积S。这意味着 2 个或更多立方体的每个体积重叠仅计算一次。该集合具体包含所有立方体的所有坐标。

我找到了几篇有关该主题的论文，介绍了完成该任务的算法。例如，本文d将问题推广到任何维度而不是琐碎的维度d=3，并将问题推广到盒子而不是立方体。这篇论文以及其他一些论文及时O(n^1.5)或稍微好一点地解决了这个问题。另一篇看起来很有前途且专门针对的3d-cubes论文解决了以下任务：O(n^4/3 log n)。但这些论文看起来相当复杂，至少对我来说是这样，我无法清楚地理解它们。

是否有任何相对简单的伪代码或过程可以遵循来实现这个想法？我正在寻找一组说明，说明如何处理立方体。任何实施也将是优秀的。甚至O(n^2)都O(n^3)很好。

---

目前，我的方法是计算总体积，即所有立方体的所有体积之和，然后计算每两个立方体的重叠，并从总体积中减去它。问题是每个这样的重叠可能（或可能不是）属于不同的立方体对，这意味着重叠可以由 5 个立方体共享。在这种方法中，重叠将被计算 10 次，而不是仅计算一次。所以我想也许包含排除原则可能会证明自己有用，但我不知道它具体是如何实施的。（天真地）计算每个重叠已经是O(n^2)可行的。有没有更好的方法来计算所有这些重叠？无论如何，我不认为这是一个有用的方法，可以继续沿着这些思路继续下去。

Answer 1

我用 C++ 实现了 Bentley 算法 (O(n^2 log n))。（我知道您想要 Java，但 C++ 是我工作中的主要斧头，考虑到我正在考虑向 Overmars 和 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";
}