heightmap_cuboid_extract.cpp

#include <iostream>
#include <vector>
#include <cassert>
#include <set>
#include <tuple>
#include <random>
#include <memory.h>
#include <limits.h>
#include <functional>

// ALIASES

using namespace std;
using u32 = uint32_t;
using u64 = uint64_t;

// CONSTANTS

constexpr u32 log2floor(u32 x) {
    return 31u - static_cast<u32>(__builtin_clz(x));
}

constexpr u32 log2ceil(u32 x) {
    u32 log2 = log2floor(x);
    bool isNotPowerOf2 = x != (1 << log2);
    return log2 + isNotPowerOf2;
}

static_assert (log2floor(1) == 0 && log2ceil(1) == 0);
static_assert (log2floor(2) == 1 && log2ceil(2) == 1);
static_assert (log2floor(3) == 1 && log2ceil(3) == 2);
static_assert (log2floor(4) == 2 && log2ceil(4) == 2);
static_assert (log2floor(7) == 2 && log2ceil(7) == 3);
static_assert (log2floor(8) == 3 && log2ceil(8) == 3);
static_assert (log2floor(15) == 3 && log2ceil(15) == 4);
static_assert (log2floor(16) == 4 && log2ceil(15) == 4);

constexpr u32 DIM_LIMIT = 100;
constexpr u32 LOG2CEIL_DIM_LIMIT = log2ceil(DIM_LIMIT);
constexpr u32 HEIGHT_LIMIT = 100;

// TYPES AND UTILITY

struct point {
    u32 x, y;
};

struct rectangle {
    u32 minX, maxX;
    u32 minY, maxY;

    u32 sizeX() const {
        return maxX - minX + 1;
    }

    u32 sizeY() const {
        return maxY - minY + 1;
    }

    u32 area() const {
        return sizeX() * sizeY();
    }
};

constexpr bool operator<(const rectangle &l, const rectangle &r) {
    return make_tuple(l.minX, l.maxX, l.minY, l.maxY) < make_tuple(r.minX, r.maxX, r.minY, r.maxY);
}

template<class T>
struct rmq2d {

    u32 h, logh;
    u32 w, logw;
    vector<vector<T>> heightmap;

    point A[LOG2CEIL_DIM_LIMIT][DIM_LIMIT][LOG2CEIL_DIM_LIMIT][DIM_LIMIT];

    rmq2d() = default;

    rmq2d(const vector<vector<T>> &heightmap) : heightmap(heightmap) {
        h = static_cast<u32>(heightmap.size());
        logh = log2floor(h) + 1;
        w = static_cast<u32>(heightmap[0].size());
        logw = log2floor(w) + 1;

        for (u32 y=0; y < h; ++y){
            for (u32 x = 0; x<w; x++) {
                A[0][y][0][x] = {x, y};
            }
            for (u32 jx=1; jx < logw; jx++){
                u32 sz = 1 << (jx - 1);
                for (u32 x=0; x+sz < w; x++) {
                    point i1 = A[0][y][jx-1][x];
                    point i2 = A[0][y][jx-1][x+sz];
                    A[0][y][jx][x] = query(i1) < query(i2) ? i1 : i2;
                }
            }
        }

        for (u32 jy=1; jy < logh; ++jy){
            u32 sz = 1 << (jy-1);
            for (u32 y=0; y+sz < h; ++y){
                for (u32 jx=0; jx < logw; ++jx){
                    for (u32 x=0; x < w; ++x){
                        point i1 = A[jy-1][y][jx][x];
                        point i2 = A[jy-1][y+sz][jx][x];
                        A[jy][y][jx][x] = query(i1) < query(i2) ? i1 : i2;
                    }
                }
            }
        }
    }

    T query(point pos) {
        return heightmap[pos.y][pos.x];
    }

    point findLowestPoint(rectangle rect) {
        u32 logLenX = log2floor(rect.sizeX());
        u32 logLenY = log2floor(rect.sizeY());

        // we must add 1 here, otherwise this can produce an underflow
        u32 maxX2 = rect.maxX + 1 - (1 << logLenX);
        u32 maxY2 = rect.maxY + 1 - (1 << logLenY);
        // alternatively clamp underflows, but this produces invalid results
        //u32 maxX2 = rect.maxX < (1 << logLenX) ? 0 : rect.maxX - (1 << logLenX);
        //u32 maxY2 = rect.maxY < (1 << logLenY) ? 0 : rect.maxY - (1 << logLenY);

        point indexes[] {
             A[logLenY][rect.minY][logLenX][rect.minX],
             A[logLenY][maxY2][logLenX][rect.minX],
             A[logLenY][rect.minY][logLenX][maxX2],
             A[logLenY][maxY2][logLenX][maxX2]
        };

        point result = indexes[0];
        for (size_t i = 1; i < 4; ++i){
            if (query(indexes[i]) < query(result)) {
                result = indexes[i];
            }
        }
        return result;
    }
};

static rmq2d<u64> rmq;
static set<rectangle> cac;
static vector<vector<u64>> heightmap(DIM_LIMIT, vector<u64>(DIM_LIMIT, 0));

static rectangle resultBase;
static u64 resultHeight = 0;
static u64 resultVolume = 0;
static size_t numOfQueries = 0;

static bool isCuboidFillable(rectangle base, u64 height) {
    for (u32 y = base.minY; y <= base.maxY; ++y) {
        for (u32 x = base.minX; x <= base.maxX; ++x) {
            if (heightmap[y][x] < height) {
                return false;
            }
        }
    }
    return true;
}

static void printIntermediateResult() {
    cout << resultBase.minX << "-" << resultBase.maxX << ", " << resultBase.minY << "-" << resultBase.maxY;
    cout << " h=" << resultHeight <<  ": ";
    cout << resultBase.sizeX() << "*" << resultBase.sizeY() << "*" << resultHeight << "=" << resultVolume << endl;
}

static void findLargestCuboid_rmq2d(rectangle rect){
    if (rect.minX >= rect.maxX || rect.minY >= rect.maxY) return;

    if (!cac.insert(rect).second)
        return;
    ++numOfQueries;

    const auto p = rmq.findLowestPoint(rect);
    u64 h = heightmap[p.y][p.x];
    u64 volume = h * rect.area();
    if (volume > resultVolume){
        resultBase = rect;
        resultHeight = h;
        resultVolume = volume;
        printIntermediateResult();
    }

    findLargestCuboid_rmq2d({p.x + 1, rect.maxX, rect.minY, rect.maxY});
    findLargestCuboid_rmq2d({rect.minX, p.x, rect.minY, rect.maxY});
    findLargestCuboid_rmq2d({rect.minX, rect.maxX, p.y + 1, rect.maxY});
    findLargestCuboid_rmq2d({rect.minX, rect.maxX, rect.minY, p.y});
}

static void findLargestCuboid_stupid(rectangle rect) {
    for (u32 h = HEIGHT_LIMIT; h-- != 0;) {
        // loop over all lower bounds
        for (u32 minY = rect.minY; minY < rect.maxY; ++minY) {
            for (u32 minX = rect.minX; minX < rect.maxX; ++minX) {
                // loop over all upper bounds
                for (u32 maxY = minY; maxY < rect.maxY; ++maxY) {
                    for (u32 maxX = minX; maxX < rect.maxX; ++maxX) {

                        const rectangle base{minX, maxX, minY, maxY};
                        const u64 volume = base.area() * h;
                        if (volume > resultVolume && isCuboidFillable(base, h)) {
                            resultBase = base;
                            resultHeight = h;
                            resultVolume = volume;
                            printIntermediateResult();
                        }

                    }
                }
            }
        }
    }
}

static void initHeightmap(size_t width, size_t height) {
    mt19937 rng{random_device()()};
    uniform_int_distribution<u32> distr{0, HEIGHT_LIMIT - 1};

    for (size_t y = 0; y < height; ++y){
        for (size_t x=0; x < width; ++x){
            heightmap[y][x] = distr(rng);
        }
    }
}

int main() {
    size_t width = heightmap[0].size();
    size_t height = heightmap.size();
    initHeightmap(width, height);

    rectangle bounds = {0, static_cast<u32>(width - 1), 0, static_cast<u32>(height - 1)};

    {
        resultVolume = 0;
        cout << "RUNNING RMQ2D ALGORITHM" << endl;
        rmq = rmq2d<u64>(heightmap);
        findLargestCuboid_rmq2d(bounds);
        cout << "numOfQueries:" << numOfQueries << endl;
        if (!isCuboidFillable(resultBase, resultHeight)) {
            cout << "INVALID RESULT!!!:" << endl;
        }
        cout << endl;
    }

    {
        resultVolume = 0;
        cout << "RUNNING STUPID ALGORITHM" << endl;
        findLargestCuboid_stupid(bounds);
        if (!isCuboidFillable(resultBase, resultHeight)) {
            cout << "INVALID RESULT!!!:" << endl;
        }
        cout << endl;
    }
}