Предложения по вычислению пересечений кратных выпуклых двумерных многоугольников

Question

Я пишу этот вопрос в поисках любого современного программного обеспечения или методов, которые могут быстро вычислить пересечение N 2D многоугольники (выпуклые оболочки проектируемых выпуклых многогранников) и M 2D полигоны, где обычно N >> M, N может быть в порядке или по крайней мере 1М полигонов и N в порядке 50к. Я искал некоторое время, но продолжаю придумывать тот же ответ, показанный ниже.

Используйте повышение и цикл, чтобы

1. вычислить проекцию многогранника (не узкого места)
2. вычислить выпуклую оболочку указанного многогранника (узкое место)
3. вычислить пересечение спроецированного многогранника и существующего 2D-многоугольника (главное узкое место).

Этот цикл повторяется NK времена, когда обычно K << M, а также K среднее число двумерных многоугольников, пересекающих один спроецированный многогранник. Это сделано, чтобы уменьшить количество вычислений.

Проблема в том, что если у меня есть N=262144 а также M=19456 это занимает около 129 секунд (при многопоточности многогранником), и это должно быть сделано около 300 раз. В идеале я хотел бы сократить время вычислений примерно до 1 секунды для указанных выше размеров, поэтому мне было интересно, может ли кто-нибудь помочь указать какое-то программное обеспечение или литературу, которые могли бы повысить эффективность.

[РЕДАКТИРОВАТЬ]

Запрос @ sehe Я выкладываю самые важные части кода. Я не скомпилировал его, так что это просто для того, чтобы понять суть... этот код предполагает, что есть воксели и пиксели, но формы могут быть любыми. Порядок точек в сетке может быть любым, но индексы расположения точек в сетке одинаковы.

#include <boost/geometry/geometry.hpp>                                
#include <boost/geometry/geometries/point.hpp>
#include <boost/geometry/geometries/ring.hpp>

const std::size_t Dimension = 2;
typedef boost::geometry::model::point<float, Dimension, boost::geometry::cs::cartesian> point_2d;
typedef boost::geometry::model::polygon<point_2d, false /* is cw */, true /* closed */> polygon_2d;
typedef boost::geometry::model::box<point_2d> box_2d;

std::vector<float> getOverlaps(std::vector<float> & projected_grid_vx, // projected voxels
                   std::vector<float> & pixel_grid_vx, // pixels
                   std::vector<int> & projected_grid_m, // number of voxels in each dimension
                   std::vector<int> & pixel_grid_m, // number of pixels in each dimension
                   std::vector<float> & pixel_grid_omega, // size of the pixel grid in cm
                   int projected_grid_size, // total number of voxels
                   int pixel_grid_size) { // total number of pixels

  std::vector<float> overlaps(projected_grid_size * pixel_grid_size);
  std::vector<float> h(pixel_grid_m.size());
  for(int d=0; d < pixel_grid_m.size(); d++) {
    h[d] = (pixel_grid_omega[2*d+1] - pixel_grid_omega[2*d]) / pixel_grid_m[d];
  }
  for(int i=0; i < projected_grid_size; i++){

    std::vector<float> point_indices(8);

    point_indices[0] = i;
    point_indices[1] = i + 1;
    point_indices[2] = i + projected_grid_m[0];
    point_indices[3] = i + projected_grid_m[0] + 1;
    point_indices[4] = i + projected_grid_m[0]       * projected_grid_m[1];
    point_indices[5] = i + projected_grid_m[0]       * projected_grid_m[1] + 1;
    point_indices[6] = i + (projected_grid_m[1] + 1) * projected_grid_m[0];
    point_indices[7] = i + (projected_grid_m[1] + 1) * projected_grid_m[0] + 1;

    std::vector<float> vx_corners(8 * projected_grid_m.size());

    for(int vn = 0; vn < 8; vn++) {
      for(int d = 0; d < projected_grid_m.size(); d++) {
    vx_corners[vn + d * 8] = projected_grid_vx[point_indices[vn] + d * projeted_grid_size];
      }
    }

    polygon_2d proj_voxel;
    for(int vn = 0; vn < 8; vn++) {
      point_2d poly_pt(vx_corners[2 * vn], vx_corners[2 * vn + 1]);
      boost::geometry::append(proj_voxel, poly_pt);
    }
    boost::geometry::correct(proj_voxel);

    polygon_2d proj_voxel_hull;
    boost::geometry::convex_hull(proj_voxel, proj_voxel_hull);

    box_2d bb_proj_vox;
    boost::geometry::envelope(proj_voxel_hull, bb_proj_vox);
    point_2d min_pt = bb_proj_vox.min_corner();
    point_2d max_pt = bb_proj_vox.max_corner();
    // then get min and max indices of intersecting bins
    std::vector<float> min_idx(projected_grid_m.size() - 1), 
      max_idx(projected_grid_m.size() - 1);

    // compute min and max indices of incidence on the pixel grid
    // this is easy assuming you have a regular grid of pixels
    min_idx[0] = std::min( (float) std::max( std::floor((min_pt.get<0>() - pixel_grid_omega[0]) / h[0] - 0.5 ), 0.), pixel_grid_m[0]-1);
    min_idx[1] = std::min( (float) std::max( std::floor((min_pt.get<1>() - pixel_grid_omega[2]) / h[1] - 0.5 ), 0.), pixel_grid_m[1]-1);
    max_idx[0] = std::min( (float) std::max( std::floor((max_pt.get<0>() - pixel_grid_omega[0]) / h[0] + 0.5 ), 0.), pixel_grid__m[0]-1);
    max_idx[1] = std::min( (float) std::max( std::floor((max_pt.get<1>() - pixel_grid_omega[2]) / h[1] + 0.5 ), 0.), pixel_grid_m[1]-1);


    // iterate only over pixels which intersect the projected voxel
    for(int iy = min_idx[1]; iy <= max_idx[1]; iy++) {
      for(int ix = min_idx[0]; ix <= max_idx[0]; ix++) {

    int idx = ix + iy * pixel_grid_size[0]; // `first' index of pixel corner point

    polygon_2d pix_poly;
    for(int pn = 0; pn < 4; pn++) {
      point_2d pix_corner_pt(
                 pixel_grid_vx[idx + pn % 2 + (pn / 2) * pixel_grid_m[0]],
                 pixel_grid_vx[idx + pn % 2 + (pn / 2) * pixel_grid_m[0] + pixel_grid_size]
                 );
      boost::geometry::append(pix_poly, pix_corner_pt);
    }
    boost::geometry::correct( pix_poly );

    //make this into a convex hull since the order of the point may be any
    polygon_2d pix_hull;
    boost::geometry::convex_hull(pix_poly, pix_hull);

    // on to perform intersection
    std::vector<polygon_2d> vox_pix_ints;
    polygon_2d vox_pix_int;

    try {
      boost::geometry::intersection(proj_voxel_hull, pix_hull, vox_pix_ints);
    } catch ( std::exception e ) {
      // skip since these may coincide at a point or line
      continue;
    }

    // both are convex so only one intersection expected
    vox_pix_int = vox_pix_ints[0];

    overlaps[i + idx * projected_grid_size] = boost::geometry::area(vox_pix_int);

      }
    } // end intersection for


  } //end projected_voxel for

  return overlaps;

}

Answer 1

Вы можете создать отношение полигона к ограничительной рамке:

Это можно сделать в вычислительном отношении один раз, чтобы получить среднее значение поли-площади к ВВ R постоянная. Или вы можете сделать это с геометрией, используя окружность, ограниченную BB, поскольку вы используете только спроецированный многогранник:

R = 0.0;
count = 0;
for (each poly) {
    count++;
    R += polyArea / itsBoundingBoxArea;
}
R = R/count;

Затем вычислите сумму пересечения ограничивающих прямоугольников.

Sbb = 0.0;
for (box1, box2 where box1.isIntersecting(box2)) {
  Sbb += box1.intersect(box2);
}

Затем:

Approximation = R * Sbb

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

В качестве альтернативы, если вы можете найти область полигонов быстрее, чем ее корпус, вы можете использовать фактическую вычисленную среднюю область поли. Это дало бы вам приличное приближение, избегая как пересечения поли, так и переноса.

Answer 2

Хм, проблема кажется похожей на "обнаружение столкновений" в игровых движках. Или "потенциально видимые множества".

Хотя я не знаю много о текущем состоянии дел, я помню, что оптимизация заключалась в том, чтобы заключать объекты в сферы, поскольку проверка перекрытий между сферами (или кругами в 2D) действительно дешева. Чтобы ускорить проверки на столкновения, объекты часто помещались в поисковые структуры (например, сферическое дерево (круговое дерево в 2D-случае)). В основном, организация пространства в иерархическую структуру, чтобы быстро выполнять запросы на перекрытия.

Поэтому мое предложение сводится к следующему: попробуйте взглянуть на алгоритмы обнаружения столкновений в игровых движках.

Answer 3

Предположение Я предполагаю, что вы имеете в виду "пересечения", а не пересечение. Кроме того, это не ожидаемый вариант использования, что большинство отдельных полисов из M а также N будет перекрываться в то же время. Если это предположение верно, то:

Ответ

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

Вот вики для quadtree:

Quadtree Wiki

Октри можно использовать в 3D.

На самом деле это даже не должно быть октри. Вы можете получить те же результаты с любым разделом пространства. Вы можете найти максимальное разделение полиса (давайте назовем это S). И создать сказать S/10 космические перегородки. Тогда у вас будет 10 отдельных пространств для параллельного выполнения. Мало того, что это будет одновременно, но больше не будет M * N время, поскольку не каждый поли должен сравниваться с любым другим поли.