建议计算多个凸2D多边形的交叉点

我正在为任何可以快速计算N 2D多边形（投影凸Polyhedrons的凸壳）的交叉的任何最先进软件或方法钓鱼。M 2D多边形通常是N >> M。N可以按顺序或至少为1M多边形，而N则以50k为单位。我现在搜索了一段时间，但是我一直想出下面显示的相同答案。

使用Boost和循环

1. 计算多面体的投影（不是瓶颈）
2. 计算上述多面体（瓶颈）的凸壳
3. 计算预测的多面体和现有2D多边形（主要瓶颈）的交点。

此循环重复NK次，其中通常为K << M，而K是与单个投影多面体相交的2D多边形的平均数量。这样做是为了减少计算数量。

问题是，如果我拥有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;
}