Computing the volume of the union of axis-aligned cubes

Question

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.

Answer 1

Here's some Python (sorry, didn't notice the Java tag) implementing user3386109's suggestion. This algorithm is O(n³ log n). We could get down to O(n³) by sorting the events for all cubes once and extracting the sorted sub-sequence that we need each time, but perhaps this is good enough.

import collections

Interval = collections.namedtuple("Interval", ["lower", "upper"])
Cube = collections.namedtuple("Cube", ["x", "y", "z"])


def length(interval):
    return interval.upper - interval.lower


def length_of_union(intervals):
    events = []
    for interval in intervals:
        events.append((interval.lower, 1))
        events.append((interval.upper, -1))
    events.sort()
    previous = None
    overlap = 0
    total = 0
    for x, delta in events:
        if overlap > 0:
            total += x - previous
        previous = x
        overlap += delta
    return total


def all_boundaries(intervals):
    boundaries = set()
    for interval in intervals:
        boundaries.add(interval.lower)
        boundaries.add(interval.upper)
    return sorted(boundaries)


def subdivide_at(interval, boundaries):
    lower = interval.lower
    for x in sorted(boundaries):  # Resort is O(n) due to Timsort.
        if x < lower:
            pass
        elif x < interval.upper:
            yield Interval(lower, x)
            lower = x
        else:
            yield Interval(lower, interval.upper)
            break


def volume_of_union(cubes):
    x_boundaries = all_boundaries(cube.x for cube in cubes)
    y_boundaries = all_boundaries(cube.y for cube in cubes)
    sub_problems = collections.defaultdict(list)
    for cube in cubes:
        for x in subdivide_at(cube.x, x_boundaries):
            for y in subdivide_at(cube.y, y_boundaries):
                sub_problems[(x, y)].append(cube.z)
    return sum(
        length(x) * length(y) * length_of_union(z_intervals)
        for ((x, y), z_intervals) in sub_problems.items()
    )


# Test code below.

import random

Point = collections.namedtuple("Point", ["x", "y", "z"])


def mean(sequence):
    n = 0
    x_bar = 0
    for x in sequence:
        n += 1
        x_bar += (x - x_bar) / n
    return x_bar


def random_interval():
    a = random.random()
    b = random.random()
    return Interval(min(a, b), max(a, b))


def in_interval(x, interval):
    return interval.lower <= x <= interval.upper


test_intervals = [random_interval() for i in range(10)]
sample_coordinates = [random.random() for i in range(1000000)]
sampled_length = mean(
    any(in_interval(x, interval) for interval in test_intervals)
    for x in sample_coordinates
)
print(length_of_union(test_intervals), sampled_length)


def random_cube():
    return Cube(random_interval(), random_interval(), random_interval())


def in_cube(point, cube):
    return (
        in_interval(point.x, cube.x)
        and in_interval(point.y, cube.y)
        and in_interval(point.z, cube.z)
    )


test_cubes = [random_cube() for i in range(10)]
sample_points = [
    Point(random.random(), random.random(), random.random()) for i in range(1000000)
]
sampled_volume = mean(
    any(in_cube(point, cube) for cube in test_cubes) for point in sample_points
)
print(volume_of_union(test_cubes), sampled_volume)

Answer 2

We can solve this problem recursively, by defining an algorithm to solve the problem for D-dimensional cubes (where D=1 are lines, 2 are squares, 3 are cubes, etc.). I assume cubes are a struct with two points (start and end) that define their boundaries.

Note: I didn't run this code, but it should work (up to silly python mistakes). Complexity is expected to be O((N log N)^D) where D is the number of dimensions, and the code should be straight forward.

# Define a recursive computation, working on a given set of cubes and
# an index of the current dimension.
def find_union_area(cubes, dim_index=0):
  if len(cubes) == 0:
    return 0

  # We want to have a list of tuples of the form `(x_val, side, cube)` where
  # side is either start or end.
  starts = [{
    'val': cube.start[dim_index],
    'side': 'start',
    'cube': cube
  } for cube in cubes]
  ends = [{
    'val': cube.start[dim_index],
    'side': 'start',
    'cube': cube
  } for cube in cubes]

  # Sort all start and end locations of the cubes along the current axis.  
  locations = sorted(starts + ends, key=lambda l: l['val'])

  # Iterate the above list, each time on the interval between x_i and x_(i+1)
  # to compute the overall volume/area minus overlap.
  result = 0

  curr = locations[0]
  assert cur['side'] == 'start'

  # Track the cubes in the current interval
  current_cubes = {loc.cube}

  for loc in locations[1:]:
    prev = curr
    curr = loc
 
    if curr['side'] == 'end':
      current_cubes.remove(curr['cube'])
    else:
      current_cubes.add(curr['cube'])

    # If the current interval has 0 length, it is non-interesting
    # This would happen if two cubes start/end at the same place along the axis
    if curr['val'] == prev['val']:
      continue

    # If the current interval is empty it is also non-interesting
    if len(current_cubes) == 0:
      continue
    
    # If we are in 1-D (i.e. this is the last dimension), the size of
    # the interval is simple
    if dim_index == len(curr.cube.start):
      result += curr['val'] - prev['val']
    # Otherwise we need to check the volume also using the next dimensions
    # for the current cubes, and then multiply by that
    else:
      curr_dim_size = curr['val'] - prev['val']
      next_dims_size = find_union_area(current_cubes, dim_index + 1)
      result += curr_dim_size * next_dims_size
 
  return result