Как вы можете найти кубоид с наибольшим объемом в карте высот? (с низкой сложностью)

Question

Мне нужно найти кубоид с наибольшим объемом, содержащийся в 2D-карте высот. Карта высот представляет собой массив размером w*d, где w - это ширина, h - это высота, а d - это глубина. В C это будет выглядеть следующим образом:

unsigned heightmap[w][d]; // all values are <= h

Я уже знаю, что существует наивный алгоритм, который может решить эту проблему со сложностью O(w*d*h). Тем не менее, я подозреваю, что существует более оптимальный метод. Это работает следующим образом, в псевдокоде Pythoni c:

resultRectangle = None
resultHeight = None
resultVolume = -1

# iterate over all heights
for loopHeight in range(0, h):
    # create a 2D bitmap from our heightmap where a 1 represents a height >= loopHeight
    bool bitmap[w][d]
    for x in range(0, w):
        for y in range(0, d):
            bitmap[x][y] = heightmap[x][y] >= loopHeight

    # obtain the greatest-volume cuboid at this particular height
    maxRectangle = maxRectangleInBitmap(bitmap)
    volume = maxRectangle.area() * loopHeight

    # compare it to our current maximum and replace it if we found a greater cuboid
    if volume > resultVolume:
        resultHeight = loopHeight
        resultVolume = volume
        resultRectangle = maxRectangle

resultCuboid = resultRectangle.withHeight(resultHeight)

Нахождение наибольшей области из всех 1 в прямоугольнике - известная проблема с O(1) сложностью на пиксель или O(w*d) в нашем дело. Таким образом, общая сложность наивного подхода составляет O(w*h*d).

Так что, как я уже говорил, мне было интересно, сможем ли мы преодолеть эту сложность. Возможно, мы сможем снизить его до O(w*d * log(h)), более разумно просматривая высоты, вместо того, чтобы «грубо насиловать» их всех.

Ответ на этот вопрос Найдите самый большой кубоид, содержащий только 1 в NxNxN двоичный массив Евгения Клюева, похоже, использует аналогичный подход, но он ошибочно (?) предполагает, что объемы, которые мы найдем на этих высотах, образуют унимодальную функцию. Если бы это было так, мы могли бы использовать Поиск Золотого сечения, чтобы выбирать высоты более разумно, но я не думаю, что мы можем.

Answer 1

ОК, еще один дубль. Большая часть «мяса» находится в функции go. Он использует ту же концепцию «расщепления», что и в моем другом ответе, но использует нисходящее динамическое программирование c с запоминанием. rmq2d реализует 2D Range Minimum Query. для размера 1000x1000 это занимает около 30 секунд (при использовании 3 ГБ памяти).

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

using namespace std;

constexpr int ilog2(int x){
    return 31 - __builtin_clz(x);
}

const int MAX_DIM = 100;




template<class T>
struct rmq2d{
    struct point{
        int x,y;

        point():x(0),y(0){}
        point(int x,int y):x(x),y(y){}
    };
    typedef point array_t[MAX_DIM][ilog2(MAX_DIM)+1][MAX_DIM];

    int h, logh;
    int w, logw;
    vector<vector<T>> v;

    array_t *A;

    rmq2d(){A=nullptr;}

    rmq2d &operator=(const rmq2d &other){
        assert(sizeof(point)==8);
        if(this == &other) return *this;
        if(!A){
            A = new array_t[ilog2(MAX_DIM)+1];
        }
        v=other.v;
        h=other.h;
        logh = other.logh;
        w=other.w;
        logw=other.logw;
        memcpy(A, other.A, (ilog2(MAX_DIM)+1)*sizeof(array_t));
        return *this;
    }

    rmq2d(const rmq2d &other){
        A = nullptr;
        *this = other;
    }

    ~rmq2d(){
        delete[] A;
    }


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

    rmq2d(vector<vector<T>> &v) : v(v){
        A = new array_t[ilog2(MAX_DIM)+1];
        h = (int)v.size();
        logh = ilog2(h) + 1;
        w = (int)v[0].size();
        logw = ilog2(w) + 1;

        for(int y=0; y<h; ++y){
            for(int x=0;x<w;x++) A[0][y][0][x] = {x, y};

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

                    }
                }
            }
        }
    }

    point pos_q(int x1, int x2, int y1, int y2){
        assert(A);
        int lenx = ilog2(x2 - x1);
        int leny = ilog2(y2 - y1);

        point idxs[] = {
                 A[leny][y1][lenx][x1],
                 A[leny][y2-(1<<leny)][lenx][x1],
                 A[leny][y1][lenx][x2-(1<<lenx)],
                 A[leny][y2-(1<<leny)][lenx][x2-(1<<lenx)]
                };
        point ret = idxs[0];
        for(int i=1; i<4; ++i){
            if(query(ret) > query(idxs[i])) ret = idxs[i];
        }
        return ret;

    }

    T val_q(int x1, int x2, int y1, int y2){
        point pos = pos_q(x1,x2,y1,y2);
        return v[pos.y][pos.x];
    }
};

rmq2d<long long> rmq;

set<tuple<int, int, int ,int>> cac;
vector<vector<long long>> v(MAX_DIM-5,vector<long long>(MAX_DIM-5,0));

long long ret = 0;
int nq = 0;

void go(int x1, int x2, int y1, int y2){
    if(x1 >= x2 || y1>=y2) return;

    if(!cac.insert(make_tuple(x1,y1,x2,y2)).second) return;
    ++nq;

    auto p = rmq.pos_q(x1, x2, y1, y2);
    long long cur = v[p.y][p.x]*(x2-x1)*(y2-y1);
    if(cur > ret){
        cout << x1 << "-" << x2 << ", " << y1 << "-" << y2 << " h=" << v[p.y][p.x] <<  " :" << cur << endl;
        ret = cur;
    }

    go(p.x+1, x2, y1, y2);
    go(x1, p.x, y1, y2);
    go(x1, x2, p.y+1, y2);
    go(x1, x2, y1, p.y);
}


int main(){
    int W = (int)v[0].size();
    int H=(int)v.size();
    for(int y=0; y<H;++y){
        for(int x=0; x<W; ++x){
            v[y][x] = rand()%10000;
        }
    }
    rmq = rmq2d<long long>(v);
    go(0,W, 0, H);
    cout << "nq:" << nq << endl;
}

Answer 2

Вот идея со значительным допущением. псевдокод:

P <- points from heightmap sorted by increasing height.
R <- set of rectangles. All maximal empty sub-rectangles for the current height.
R.add(Rectangle(0,0,W,H)
result = last_point_in(P).height()
foreach(p in P):
   RR <- rectangles from R that overlap P (can be found in O(size(RR)), possibly with some logarithmic factors)
   R = R - RR
   foreach(r in RR)
       result = max(result, r.area() * p.height())
       split up r, adding O(1) new rectangles to R.
return result

Предположение, которое у меня есть внутреннее чувство, но я не могу доказать, состоит в том, что RR будет в среднем иметь размер O (1).

Править : чтобы уточнить «расщепление», если мы разбиваем в точке p:

AAAAADFFF
AAAAADFFF
AAAAADFFF
BBBBBpGGG
CCCCCEHHH
CCCCCEHHH

Мы создаем новые прямоугольники, состоящие из: AB C, CEH, FGH, ADF, и добавляем их в R.