关于计算多个凸2D多边形的交集的建议

我正在为任何可以快速计算N2D多边形（投影凸多面体的凸包）和M2D多边形（通常为）的交集的最新软件或方法编写此问题N >> M。N可以是大约1M的多边形顺序，也可以是N50k的顺序。我已经搜索了一段时间，但是我一直想出下面显示的相同答案。

使用升压和循环来

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;

}