rectangle-union-area/main.cpp

#include <vector>
#include <assert.h>

#include <benchmark/benchmark.h>

#include "geometry.h"

struct Rectangle {
    Point bottom_left;
    Point top_right;

    // only for algorithm
    Index i_left;
    Index i_right;
};

struct RectCoord {
    Coordinate coord;
    Rectangle* rect; // guaranteed to be valid at all times
    bool operator<(const RectCoord &other) {
        return coord < other.coord;
    }
};

Unit get_area_union(std::vector<Rectangle> rectangles) {
    // do some input sanity checks
    for(const auto &rect : rectangles) {
        assert(rect.bottom_left.x <= rect.top_right.x);
        assert(rect.bottom_left.y <= rect.top_right.y);
    }
    // entry i will describe coverage of the interval (x_points[i], x_points[i+1])
    std::vector<unsigned> coverage_numbers;
    coverage_numbers.reserve(2*rectangles.size() -1);

    //sorted vector of all x coordinates of all rectangles
    std::vector<RectCoord> x_points;
    x_points.reserve(2* rectangles.size());
    for(auto &rect : rectangles) {
        x_points.push_back({rect.bottom_left.x, &rect});
        x_points.push_back({rect.top_right.x, &rect});
    }
    std::sort(x_points.begin(), x_points.end());

    // preprocessing for constant rectangle lookup during algorithm
    // technically not necessary to achieve quadratic running time
    for(Index i = 0 ; i < x_points.size() ; ++i) {
        if(x_points[i].rect->bottom_left.x == x_points[i].coord) {
            x_points[i].rect->i_left = i;
        } else {
            assert(x_points[i].rect->top_right.x == x_points[i].coord);
            x_points[i].rect->i_right = i;
        }
    }

    // prepare vector of y-coordinates
    std::vector<RectCoord> y_points;
    y_points.reserve(2*rectangles.size());
    for(auto &rect : rectangles) {
        y_points.push_back({rect.bottom_left.y, &rect});
        y_points.push_back({rect.top_right.y, &rect});
    }
    std::sort(y_points.begin(), y_points.end());

    // run sweepline
    Unit current_cross_section_length = 0;
    Unit total_area = 0;
    for(Index y_index = 0 ; y_index < y_points.size() ; ++y_index) {
        // first, update cross section
        if (y_points[y_index].rect->bottom_left.y == y_points[y_index].coord) {
            // rectangle is starting now, add its cross section
            for (Index index = y_points[y_index].rect->i_left; index < y_points[y_index].rect->i_right; ++index) {
                ++coverage_numbers[index];
                if (coverage_numbers[index] == 1) {
                    current_cross_section_length += (x_points[index + 1].coord - x_points[index].coord);
                }
            }
        } else {
            assert(y_points[y_index].rect->top_right.y == y_points[y_index].coord);
            // rectangle stopping, remove its cross section
            for (Index index = y_points[y_index].rect->i_left; index < y_points[y_index].rect->i_right; ++index) {
                --coverage_numbers[index];
                if (coverage_numbers[index] == 0) {
                    current_cross_section_length -= (x_points[index + 1].coord - x_points[index].coord);
                }
            }
        }
        //cross section is now up to date
        if(y_index + 1 < y_points.size()) {
            // add area up to next y point
            total_area += current_cross_section_length * (y_points[y_index+1].coord - y_points[y_index].coord);
        }
    }
    // sanity check that nothing is covered when arriving at the end
    for(auto coverage : coverage_numbers) {
        assert(coverage == 0);
    }
    return total_area;
}

// some benchmark tests using github.com/google/benchmark

std::vector<Rectangle> get_random_instance(unsigned num_rects, Coordinate range) {
    std::vector<Rectangle> rects;
    rects.reserve(num_rects);
    for(unsigned i = 0 ; i < num_rects; ++i) {
        // technically, this might not be (perfectly) evenly distributed, but it will suffice for our purposes anyways
        Coordinate x1 = rand() % range;
        Coordinate x2 = rand() % range;
        Coordinate y1 = rand() % range;
        Coordinate y2 = rand() % range;
        if (x2 < x1) { int tmp = x1; x1 = x2; x2 = tmp; };
        if (y2 < y1) { int tmp = y1; y1 = y2; y2 = tmp; };
        rects.push_back(
                {{x1,y1},{x2, y2}, 0,0}
        );
    }
    return rects;
}

static void BM_area_computation(benchmark::State& state) {
    std::vector<Rectangle> rects;
    for (auto _ : state) {
        state.PauseTiming();
        rects = get_random_instance(state.range(0), state.range(1));
        state.ResumeTiming();
        get_area_union(rects);
    }
    state.SetComplexityN(state.range(0));
}

// add benchmark
BENCHMARK(BM_area_computation)->ArgsProduct({benchmark::CreateRange(1,1<<16,4),{1<<10, 1<<15, 1<<20}})->Complexity(benchmark::oNSquared);

// run benchmark as main()
BENCHMARK_MAIN();