#include #include #include #include "geometry.h" Unit get_area_union(std::vector 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 coverage_numbers; coverage_numbers.reserve(2*rectangles.size() -1); //sorted vector of all x coordinates of all rectangles std::vector 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 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 get_random_instance(unsigned num_rects, Coordinate range) { std::vector 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 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();