Максимальный объем ловушки дождевой воды в 3D

Классический алгоритм в 2D-версии обычно описывается как

Учитывая n неотрицательных целых чисел, представляющих карту высот, где ширина каждого бара равна 1, вычислите, сколько воды оно может поймать после дождя.

Например, учитывая ввод

[0,1,0,2,1,0,1,3,2,1,2,1] 

возвращаемое значение будет

6

enter image description here

Алгоритм, который я использовал для решения вышеупомянутой двумерной задачи, -

int trapWaterVolume2D(vector<int> A) {
    int n = A.size();
    vector<int> leftmost(n, 0), rightmost(n, 0);

    //left exclusive scan, O(n), the highest bar to the left each point
    int leftMaxSoFar = 0;
    for (int i = 0; i < n; i++){
        leftmost[i] = leftMaxSoFar;
        if (A[i] > leftMaxSoFar) leftMaxSoFar = A[i];
    }


    //right exclusive scan, O(n), the highest bar to the right each point
    int rightMaxSoFar = 0;
    for (int i = n - 1; i >= 0; i--){
        rightmost[i] = rightMaxSoFar;
        if (A[i] > rightMaxSoFar) rightMaxSoFar = A[i];
    }

    // Summation, O(n)
    int vol = 0;
    for (int i = 0; i < n; i++){
        vol += max(0, min(leftmost[i], rightmost[i]) - A[i]);
    }
    return vol;
}

Мой вопрос заключается в том, как сделать вышеуказанный алгоритм доступным для 3D-версии проблемы, чтобы вычислить максимум воды, захваченной в реальном 3D-ландшафте. т.е. реализовать

int trapWaterVolume3D(vector<vector<int> > A);

Пример графика:

enter image description here

Мы знаем высоту в каждой точке (x, y), и цель состоит в том, чтобы вычислить максимальный объем воды, который может быть захвачен в форме. Любые мысли и ссылки приветствуются.

Ответы

Ответ 1

Для каждой точки на местности рассмотрите все пути от этой точки до границы местности. Уровень воды будет минимальным максимальными высотами точек этих путей. Чтобы найти это, нам нужно выполнить слегка модифицированный алгоритм Дейкстры, заполняя матрицу уровня воды, начиная с границы.

For every point on the border set the water level to the point height
For every point not on the border set the water level to infinity
Put every point on the border into the set of active points
While the set of active points is not empty:
    Select the active point P with minimum level
    Remove P from the set of active points
    For every point Q adjacent to P:
        Level(Q) = max(Height(Q), min(Level(Q), Level(P)))
        If Level(Q) was changed:
            Add Q to the set of active points

Ответ 2

user3290797 "слегка модифицированный алгоритм Дейкстры" ближе к алгоритму Прима, чем к Dijkstra. В минимальных связующих частях дерева мы создаем граф с одной вершиной на одну плиту, одну вершину для внешней стороны и ребрами с весами, равными максимальной высоте их двух смежных плит (внешняя сторона имеет высоту "минус бесконечность" ).

Учитывая путь в этом графе во внешнюю вершину, максимальный вес ребра в пути - это высота, которую должна достичь вода, чтобы вырваться по этому пути. Соответствующее свойство минимального связующего дерева состоит в том, что для каждой пары вершин максимальный вес ребра в пути в охватывающем дереве является минимальным возможным среди всех путей между этими вершинами. Таким образом, минимальное остовное дерево описывает наиболее экономичные пути выхода воды, и высоты воды могут быть извлечены в линейном времени с одним обходом.

В качестве бонуса, поскольку график является планарным, существует алгоритм линейного времени для вычисления минимального остовного дерева, состоящий из чередующихся проходов и упрощений Boruvka. Это улучшает время работы O (n log n) Prim.

Ответ 3

Эта проблема очень близка к построению морфологического водораздела изображения в оттенках серого (http://en.wikipedia.org/wiki/Watershed_(image_processing)).

Один подход заключается в следующем (процесс наводнения):

  • сортировать все пиксели, увеличивая высоту.

  • работают постепенно, увеличивая возвышения, назначая метки пикселям на бассейн водосбора.

  • для нового уровня высоты вам нужно пометить новый набор пикселей:

    • Некоторые из них не помечены соседний, они образуют локальную минимальную конфигурацию и начинают новый бассейн.
    • У некоторых есть только соседи с одинаковой меткой, их можно пометить аналогично (они расширяют бассейн водосбора).
    • У некоторых есть соседи с разными метками. Они не относятся к определенному бассейну водосбора, и они определяют линии водораздела.

Вам нужно будет улучшить стандартный алгоритм водораздела, чтобы иметь возможность вычислить объем воды. Вы можете это сделать, определив максимальный уровень воды в каждом бассейне и выведите высоту земли на каждый пиксель. Уровень воды в бассейне определяется возвышением самого низкого водораздельного пикселя вокруг него.

Вы можете действовать каждый раз, когда обнаруживаете водораздельный пиксель: если соседнему бассейну еще не присвоен уровень, этот бассейн может выдерживать текущий уровень без утечки.

Ответ 4

Эта проблема может быть решена с использованием алгоритма Priority-Flood. Он был обнаружен и опубликован несколько раз за последние несколько десятилетий (и опять же другими людьми, отвечающими на этот вопрос), хотя конкретный вариант, который вы ищете, не является, насколько мне известно, в литературе.

Вы можете найти обзорную статью алгоритма и его варианты здесь. После публикации этой статьи был обнаружен еще более быстрый вариант (ссылка), а также методы для выполнения этого расчета на наборах данных триллионов ячеек (ссылка). Обсуждается метод выборочного нарушения низких/узких делений здесь. Свяжитесь со мной, если вы хотите получить копии любой из этих документов.

У меня есть репозиторий здесь со многими из вышеперечисленных вариантов; дополнительные реализации можно найти здесь.

Простой script для вычисления объема с использованием библиотеки RichDEM выглядит следующим образом:

#include "richdem/common/version.hpp"
#include "richdem/common/router.hpp"
#include "richdem/depressions/Lindsay2016.hpp"
#include "richdem/common/Array2D.hpp"

/**
  @brief  Calculates the volume of depressions in a DEM
  @author Richard Barnes ([email protected])

    Priority-Flood starts on the edges of the DEM and then works its way inwards
    using a priority queue to determine the lowest cell which has a path to the
    edge. The neighbours of this cell are added to the priority queue if they
    are higher. If they are lower, then they are members of a depression and the
    elevation of the flooding minus the elevation of the DEM times the cell area
    is the flooded volume of the cell. The cell is flooded, total volume
    tracked, and the neighbors are then added to a "depressions" queue which is
    used to flood depressions. Cells which are higher than a depression being
    filled are added to the priority queue. In this way, depressions are filled
    without incurring the expense of the priority queue.

  @param[in,out]  &elevations   A grid of cell elevations

  @pre
    1. **elevations** contains the elevations of every cell or a value _NoData_
       for cells not part of the DEM. Note that the _NoData_ value is assumed to
       be a negative number less than any actual data value.

  @return
    Returns the total volume of the flooded depressions.

  @correctness
    The correctness of this command is determined by inspection. (TODO)
*/
template <class elev_t>
double improved_priority_flood_volume(const Array2D<elev_t> &elevations){
  GridCellZ_pq<elev_t> open;
  std::queue<GridCellZ<elev_t> > pit;
  uint64_t processed_cells = 0;
  uint64_t pitc            = 0;
  ProgressBar progress;

  std::cerr<<"\nPriority-Flood (Improved) Volume"<<std::endl;
  std::cerr<<"\nC Barnes, R., Lehman, C., Mulla, D., 2014. Priority-flood: An optimal depression-filling and watershed-labeling algorithm for digital elevation models. Computers & Geosciences 62, 117–127. doi:10.1016/j.cageo.2013.04.024"<<std::endl;

  std::cerr<<"p Setting up boolean flood array matrix..."<<std::endl;
  //Used to keep track of which cells have already been considered
  Array2D<int8_t> closed(elevations.width(),elevations.height(),false);

  std::cerr<<"The priority queue will require approximately "
           <<(elevations.width()*2+elevations.height()*2)*((long)sizeof(GridCellZ<elev_t>))/1024/1024
           <<"MB of RAM."
           <<std::endl;

  std::cerr<<"p Adding cells to the priority queue..."<<std::endl;

  //Add all cells on the edge of the DEM to the priority queue
  for(int x=0;x<elevations.width();x++){
    open.emplace(x,0,elevations(x,0) );
    open.emplace(x,elevations.height()-1,elevations(x,elevations.height()-1) );
    closed(x,0)=true;
    closed(x,elevations.height()-1)=true;
  }
  for(int y=1;y<elevations.height()-1;y++){
    open.emplace(0,y,elevations(0,y)  );
    open.emplace(elevations.width()-1,y,elevations(elevations.width()-1,y) );
    closed(0,y)=true;
    closed(elevations.width()-1,y)=true;
  }

  double volume = 0;

  std::cerr<<"p Performing the improved Priority-Flood..."<<std::endl;
  progress.start( elevations.size() );
  while(open.size()>0 || pit.size()>0){
    GridCellZ<elev_t> c;
    if(pit.size()>0){
      c=pit.front();
      pit.pop();
    } else {
      c=open.top();
      open.pop();
    }
    processed_cells++;

    for(int n=1;n<=8;n++){
      int nx=c.x+dx[n];
      int ny=c.y+dy[n];
      if(!elevations.inGrid(nx,ny)) continue;
      if(closed(nx,ny))
        continue;

      closed(nx,ny)=true;
      if(elevations(nx,ny)<=c.z){
        if(elevations(nx,ny)<c.z){
          ++pitc;
          volume += (c.z-elevations(nx,ny))*std::abs(elevations.getCellArea());
        }
        pit.emplace(nx,ny,c.z);
      } else
        open.emplace(nx,ny,elevations(nx,ny));
    }
    progress.update(processed_cells);
  }
  std::cerr<<"t Succeeded in "<<std::fixed<<std::setprecision(1)<<progress.stop()<<" s"<<std::endl;
  std::cerr<<"m Cells processed = "<<processed_cells<<std::endl;
  std::cerr<<"m Cells in pits = "  <<pitc           <<std::endl;

  return volume;
}

template<class T>
int PerformAlgorithm(std::string analysis, Array2D<T> elevations){
  elevations.loadData();

  std::cout<<"Volume: "<<improved_priority_flood_volume(elevations)<<std::endl;

  return 0;
}

int main(int argc, char **argv){
  std::string analysis = PrintRichdemHeader(argc,argv);

  if(argc!=2){
    std::cerr<<argv[0]<<" <Input>"<<std::endl;
    return -1;
  }

  return PerformAlgorithm(argv[1],analysis);
}

Это должно быть прямолинейно адаптировать это к любому формату 2d массива, который вы используете

В псевдокоде следующее эквивалентно сказанному:

Let PQ be a priority-queue which always pops the cell of lowest elevation
Let Closed be a boolean array initially set to False
Let Volume = 0
Add all the border cells to PQ.
For each border cell, set the cell entry in Closed to True.
While PQ is not empty:
  Select the top cell from PQ, call it C.
  Pop the top cell from PQ.
  For each neighbor N of C:
    If Closed(N):
      Continue
    If Elevation(N)<Elevation(C):
      Volume += (Elevation(C)-Elevation(N))*Area
      Add N to PQ, but with Elevation(C)
    Else:
      Add N to PQ with Elevation(N)
    Set Closed(N)=True

Ответ 5

Чтобы выполнить проблему водопроводной воды в 3D, то есть, чтобы вычислить максимальный объем захваченной дождевой воды, вы можете сделать что-то вроде этого:

#include<bits/stdc++.h>
using namespace std;

#define MAX 10

int new2d[MAX][MAX];
int dp[MAX][MAX],visited[MAX][MAX];


int dx[] = {1,0,-1,0};
int dy[] = {0,-1,0,1};


int boundedBy(int i,int j,int k,int in11,int in22)
{
    if(i<0 || j<0 || i>=in11 || j>=in22)
        return 0;

    if(new2d[i][j]>k)
        return new2d[i][j];

    if(visited[i][j])   return INT_MAX;

    visited[i][j] = 1;

    int r = INT_MAX;
    for(int dir = 0 ; dir<4 ; dir++)
    {
        int nx = i + dx[dir];
        int ny = j + dy[dir];
        r = min(r,boundedBy(nx,ny,k,in11,in22));
    }
    return r;
}

void mark(int i,int j,int k,int in1,int in2)
{
    if(i<0 || j<0 || i>=in1 || j>=in2)
        return;

    if(new2d[i][j]>=k)
        return;

    if(visited[i][j])   return ;

    visited[i][j] = 1;

    for(int dir = 0;dir<4;dir++)
    {
        int nx = i + dx[dir];
        int ny = j + dy[dir];
        mark(nx,ny,k,in1,in2);
    }
    dp[i][j] = max(dp[i][j],k);
}

struct node
{
    int i,j,key;
    node(int x,int y,int k)
    {
        i = x;
        j = y;
        key = k;
    }
};

bool compare(node a,node b)
{
    return a.key>b.key;
}
vector<node> store;

int getData(int input1, int input2, int input3[])
 {

    int row=input1;
    int col=input2;
    int temp=0;
    int count=0;

        for(int i=0;i<row;i++)
        {
            for(int j=0;j<col;j++)
            {
                if(count==(col*row)) 
                break;
            new2d[i][j]=input3[count];
            count++;
            }
        }

    store.clear();
    for(int i = 0;i<input1;i++)
        {
            for(int j = 0;j<input2;j++)
            {
                store.push_back(node(i,j,new2d[i][j]));
            }
        }
     memset(dp,0,sizeof(dp));

        sort(store.begin(),store.end(),compare);

       for(int i = 0;i<store.size();i++)
        {
            memset(visited,0,sizeof(visited));

            int aux = boundedBy(store[i].i,store[i].j,store[i].key,input1,input2);
            if(aux>store[i].key)
            {

                 memset(visited,0,sizeof(visited));
                 mark(store[i].i,store[i].j,aux,input1,input2);
            }

        }

        long long result =0 ;

        for(int i = 0;i<input1;i++)
        {

            for(int j = 0;j<input2;j++)
            {
                result = result + max(0,dp[i][j]-new2d[i][j]);
            }
        }

      return result;

 }


int main()
{
    cin.sync_with_stdio(false);
    cout.sync_with_stdio(false);
       int n,m;
        cin>>n>>m;
       int inp3[n*m];
        store.clear();

            for(int j = 0;j<n*m;j++)
            {
                cin>>inp3[j];

            }

    int k = getData(n,m,inp3);
       cout<<k;

    return 0;
}

Ответ 6

Вот простой код для этого же -

#include<iostream>
using namespace std;
int main()
{
   int n,count=0,a[100];
   cin>>n;
   for(int i=0;i<n;i++)
   {
       cin>>a[i];
   }
   for(int i=1;i<n-1;i++)
   {
       ///computing left most largest and Right most largest element of array;
       int leftmax=0;
       int rightmax=0;
       ///left most largest
       for(int j=i-1;j>=1;j--)
       {
           if(a[j]>leftmax)
           {
              leftmax=a[j];
           }
       }
       ///rightmost largest

       for(int k=i+1;k<=n-1;k++)
       {
           if(a[k]>rightmax)
           {
              rightmax=a[k];
           }
       }
       ///computing hight of the water contained-

       int x=(min(rightmax,leftmax)-a[i]);
       if(x>0)
       {
           count=count+x;
       }

   }
   cout<<count;
   return 0;
}