algorithm - El volumen máximo de agua de lluvia atrapada en 3D

graphics computational-geometry

Question

Una pregunta de algoritmo clásico en versión 2D se describe típicamente como

Dados n enteros no negativos que representan un mapa de elevación donde el ancho de cada barra es 1, calcule cuánta agua puede atrapar después de llover.

Por ejemplo, dada la entrada

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

el valor de retorno seria

6

El algoritmo que utilicé para resolver el problema 2D anterior es

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;
}

Mi pregunta es cómo hacer que el algoritmo anterior sea extensible a la versión 3D del problema, para calcular el máximo de agua atrapada en un terreno 3D del mundo real. es decir, para implementar

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

Gráfico de muestra:

Sabemos la elevación en cada punto (x, y) y el objetivo es calcular el volumen máximo de agua que puede quedar atrapado en la forma. Cualquier pensamiento y referencias son bienvenidos.

Answer 1

Aquí está el código simple para el mismo

#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;
}

Answer 2

El "algoritmo Dijkstra ligeramente modificado" de user3290797 está más cerca del algoritmo de Prim que el de Dijkstra. En términos de árbol de expansión mínimo, preparamos una gráfica con un vértice por mosaico, un vértice para el exterior y bordes con pesos iguales a la altura máxima de sus dos mosaicos contiguos (el exterior tiene una altura "menos infinito").

Dado un camino en este gráfico hacia el vértice exterior, el peso máximo de un borde en el camino es la altura que el agua debe alcanzar para escapar por ese camino. La propiedad relevante de un árbol de expansión mínimo es que, para cada par de vértices, el peso máximo de un borde en la ruta en el árbol de expansión es el mínimo posible entre todas las rutas entre esos vértices. El árbol de expansión mínimo, por lo tanto, describe las rutas de escape más económicas para el agua, y las alturas de agua se pueden extraer en tiempo lineal con un recorrido.

Como beneficio adicional, dado que el gráfico es plano, hay un algoritmo de tiempo lineal para calcular el árbol de expansión mínimo, que consiste en alternar pases Boruvka y simplificaciones. Esto mejora el tiempo de ejecución O (n log n) de Prim.

Answer 3

Este problema está muy cerca de la construcción de la cuenca morfológica de una imagen en escala de grises ( http://en.wikipedia.org/wiki/Watershed_(image_processing) ).

Un enfoque es el siguiente (proceso de inundación):

• ordena todos los píxeles aumentando la elevación.

• trabaje de forma incremental, aumentando las elevaciones, asignando etiquetas a los píxeles por cuenca de captación.

• para un nuevo nivel de elevación, necesita etiquetar un nuevo conjunto de píxeles:

  • Algunos no tienen vecinos etiquetados, forman una configuración mínima local y comienzan una nueva cuenca de captación.
  • Algunos solo tienen vecinos con la misma etiqueta, se pueden etiquetar de manera similar (extienden una cuenca de captación).
  • Algunos tienen vecinos con diferentes etiquetas. No pertenecen a una cuenca de captación específica y definen las líneas de cuencas hidrográficas.

Deberá mejorar el algoritmo de cuenca hidrográfica estándar para poder calcular el volumen de agua. Puede hacerlo determinando el nivel máximo de agua en cada cuenca y deduciendo la altura del suelo en cada píxel. El nivel de agua en una cuenca viene dado por la elevación del píxel de cuenca más bajo que lo rodea.

Puede actuar cada vez que descubra un píxel de cuenca hidrográfica: si a una cuenca vecina todavía no se le ha asignado un nivel, esa cuenca puede soportar el nivel actual sin fugas.

Answer 4

Este problema se puede resolver utilizando el algoritmo de inundación de prioridad. Se ha descubierto y publicado varias veces en las últimas décadas (y otra vez por otras personas que responden a esta pregunta), aunque la variante específica que está buscando no está, según mi conocimiento, en la literatura.

Puede encontrar un artículo de revisión del algoritmo y sus variantes here . Desde que se publicó ese artículo, se ha descubierto una variante aún más rápida ( link ), así como métodos para realizar este cálculo en conjuntos de datos de billones de células ( link ). Un método para romper selectivamente divisiones bajas / estrechas se discute here . Póngase en contacto conmigo si desea copias de cualquiera de estos documentos.

Tengo un repositorio here con muchas de las variantes anteriores; Implementaciones adicionales se pueden encontrar here .

Una secuencia de comandos simple para calcular el volumen utilizando la biblioteca RichDEM es la siguiente:

#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);
}

Debe ser sencillo adaptar esto al formato de matriz 2d que esté utilizando

En pseudocódigo, lo siguiente es equivalente a lo anterior:

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''s 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

Answer 5

Para cada punto en el terreno, considere todos los caminos desde ese punto hasta el borde del terreno. El nivel de agua sería el mínimo de las alturas máximas de los puntos de esos caminos. Para encontrarlo, necesitamos realizar un algoritmo de Dijkstra ligeramente modificado, rellenando la matriz del nivel del agua desde el borde.

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

Answer 6

Para resolver el problema del agua del grifo en 3D, es decir, para calcular el volumen máximo de agua de lluvia atrapada, puede hacer algo como esto:

#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;
}

Answer 7

    class Solution(object):
def trapRainWater(self, heightMap):
    """
    :type heightMap: List[List[int]]
    :rtype: int
    """
    m = len(heightMap)
    if m == 0:
        return 0
    n = len(heightMap[0])
    if n == 0:
        return 0
    visited = [[False for i in range(n)] for j in range(m)]
    from Queue import PriorityQueue
    q = PriorityQueue()
    for i in range(m):
        visited[i][0] = True
        q.put([heightMap[i][0],i,0])
        visited[i][n-1] = True
        q.put([heightMap[i][n-1],i,n-1])
    for j in range(1, n-1):
        visited[0][j] = True
        q.put([heightMap[0][j],0,j])
        visited[m-1][j] = True
        q.put([heightMap[m-1][j],m-1,j])
    S = 0
    while not q.empty():
        cell = q.get()
        for (i, j) in [(1,0), (-1,0), (0,1), (0,-1)]:
            x = cell[1] + i
            y = cell[2] + j
            if x in range(m) and y in range(n) and not visited[x][y]:
                S += max(0, cell[0] - heightMap[x][y])   # how much water at the cell
                q.put([max(heightMap[x][y],cell[0]),x,y])
                visited[x][y] = True
    return S