2

我实际上在优化我的算法时遇到了一些麻烦：

我有一个磁盘（以 0 为中心，半径为 1）填充了三角形（不一定具有相同的面积/长度）。可能有大量的三角形（比如说从1k到300k三角形）

我的目标是尽快找到一个点属于哪个三角形。该操作必须重复很长时间（大约10ktimes）。

现在我使用的算法是：我正在计算每个三角形中点的重心坐标。如果第一个系数介于 0 和 1 之间，我继续。如果不是，我就停下来。然后我用同样的想法计算第二个系数，第三个系数，我对每个三角形都这样做。

我想不出一种方法来使用我正在处理光盘的事实（以及我有一个欧几里得距离来帮助我直接“瞄准”好三角形的事实）：如果我尝试计算距离我指向三角形的每个“中心”：

1）当我用重心坐标强制它时，它已经比我正在做的操作更多

2）我必须订购一个包含所有三角形的欧几里得距离的向量到我的点。

3）我绝对不能保证离我的点最近的三角形是好的三角形。

我觉得我错过了一些东西，我可以预先计算，一些东西可以帮助我在开始蛮力部分之前发现好的“区域”。

该算法已经并行化（使用 OpenMP）：我正在并行调用以下函数。

bool Triangle2D::is_in_triangle(Vector2d Point) {
    double denominator = ((Tri2D(1, 1) - Tri2D(2, 1))*(Tri2D(0, 0) - Tri2D(2, 0)) + (Tri2D(2, 0) - Tri2D(1, 0))*(Tri2D(0, 1) - Tri2D(2, 1)));
    // Computing the first coefficient
    double a = ((Tri2D(1, 1) - Tri2D(2, 1))*(Point(0) - Tri2D(2, 0)) + (Tri2D(2, 0) - Tri2D(1, 0))*(Point(1) - Tri2D(2, 1))) / denominator;
    if (a < 0 || a>1) {
        return(false);
    }
    // Computing the second coefficient
    double b = ((Tri2D(2, 1) - Tri2D(0, 1))*(Point(0) - Tri2D(2, 0)) + (Tri2D(0, 0) - Tri2D(2, 0))*(Point(1) - Tri2D(2, 1))) / denominator;
    if (b < 0 || b>1) {
        return(false);
    }
    // Computing the third coefficient
    double c = 1 - a - b;
    if (c < 0 || c>1) {
        return(false);
    }
    return(true);
}

下一步可能是查看 GPU 并行化，但我需要确保代码背后的想法足够好。

现在它大约2min30需要75k三角形和10k点，但这还不够快。

编辑：
Triangle2D使用特征矩阵来存储坐标

4

4 回答 4

2

所有留着长胡须的HPC 专业人士，请允许在这里进行一些学术上精心设计的方法，对于我们的社区成员来说，这可能会变得有趣，如果不喜欢的话，他们会觉得自己比你的专业感觉更初级您自己以及可能对更深入地了解以性能为动力的代码设计、性能调整和其他并行代码风险和好处感兴趣的人，您对自己的核心 HPC 计算经验非常了解。谢谢你。

a) 算法（原样）可以获得约 2 倍的加速 一个唾手可得的果实 + 更多惊喜尚未到来

b) 由于2geometry，其他算法可能会获得~40~80X 的加速提升

c) 最佳并行代码 + 终极性能的提示

目标：在具有 i5 8500、3GHz 、6 核、NVIDIA Quadro P400 的计算机上，三角形10k点的目标运行时间为2-3分钟（必须尝试 GPU 计算，甚至不确定是否值得）300k

在此处输入图像描述

虽然这似乎是一个漫长的旅程，但这个问题很好，值得仔细研究一下，所以请耐心等待我的最佳性能动机思维流程。

a) 算法（原样）分析：

原样使用重心坐标系是一个不错的技巧，在最佳情况下，直接实现的成本略高于（20 FLOPs + 16 MEM/REG-I/O-ops），略高于（30 FLOPs + 30 MEM/REG-I/O 操作）。

有一些改进，可以通过避免一些昂贵甚至不重要的操作来降低这些执行成本：

--------------------------------------------------------------------------------------
double denominator = ( ( Tri2D( 1, 1 )
                       - Tri2D( 2, 1 ) // -------------------------- 2x MEM + OP-1.SUB
                         ) * ( Tri2D( 0, 0 ) //---------------------        + OP-3.MUL
                             - Tri2D( 2, 0 ) //--------------------- 2x MEM + OP-2.SUB
                               ) + ( Tri2D( 2, 0 ) //---------------        + OP-7.ADD
                                   - Tri2D( 1, 0 ) //--------------- 2x MEM + OP-4.SUB
                                     ) * ( Tri2D( 0, 1 ) //---------        + OP-6.MUL
                                         - Tri2D( 2, 1 ) //--------- 2x MEM + OP-5.SUB
                                           )
                       );
// Computing the first coefficient ------------------------------------------------------------------------------------------------------
double           a = ( ( Tri2D( 1, 1 )
                       - Tri2D( 2, 1 )  //-------------------------- 2x MEM + OP-8.SUB
                         ) * ( Point(0)   //------------------------        + OP-A.MUL
                             - Tri2D( 2, 0 ) //--------------------- 2x MEM + OP-9.SUB
                               ) + ( Tri2D( 2, 0 ) //---------------        + OP-E.ADD
                                   - Tri2D( 1, 0 ) //--------------- 2x MEM + OP-B.SUB
                                     ) * ( Point(1) //--------------        + OP-D.MUL
                                         - Tri2D( 2, 1 ) //--------- 2x MEM + OP-C.MUL
                                           )
                       ) / denominator; //-------------------------- 1x REG + OP-F.DIV //----------- MAY DEFER THE MOST EXPENSIVE DIVISION UNTIL a third coeff is indeed first needed, if ever------------[3]
if (a < 0 || a>1) { // ----------------------------------------------------------------------------- a < 0 ~~    ( sign( a ) * sign( denominator ) ) < 0
    return(false); // ------------------------------------------------------------------------------ a > 1 ~~  ||        a >         denominator
}
// Computing the second coefficient
double           b = ( ( Tri2D( 2, 1 ) - Tri2D( 0, 1 ) ) //--------- 2x MEM + OP-16.SUB
                     * (      Point(0) - Tri2D( 2, 0 ) ) //--------- 2x MEM + OP-17.SUB + OP-18.MUL
                     + ( Tri2D( 0, 0 ) - Tri2D( 2, 0 ) ) //--------- 2x MEM + OP-19.SUB + OP-22.ADD
                     * (      Point(1) - Tri2D( 2, 1 ) ) //--------- 2x MEM + OP-20.SUB + OP-21.MUL
                       ) / denominator; //-------------------------- 1x REG + OP-23.DIV //---------- MAY DEFER THE MOST EXPENSIVE DIVISION UNTIL a third coeff is indeed first needed, if ever -----------[3]
if (b < 0 || b>1) { // ----------------------------------------------------------------------------- b < 0 ~~   ( sign( b ) * sign( denominator ) ) < 0
    return(false); // ------------------------------------------------------------------------------ b > 1 ~~ ||        b >         denominator
}
// Computing the third coefficient
double c = 1 - a - b; // ------------------------------------------- 2x REG + OP-24.SUB + OP-25.SUB
//         1 -(a - b)/denominator; //--------------------------------------------------------------- MAY DEFER THE MOST EXPENSIVE DIVISION EXECUTED BUT HERE, IFF INDEED FIRST NEEDED <---HERE <----------[3]

  • 重复的重新评估（出现在原始版本中）可能会通过手动分配/重用明确地制作出来，但是，一个好的优化编译器可能会在使用-O3强制优化标志的情况下将它们驱逐出去。

  • 即使是这个最容易实现的果实，也要毫不犹豫地对最昂贵的部分进行打磨。

在此处输入图像描述

//------------------------------------------------------------------
double Tri2D_11_sub_21 = ( Tri2D( 1, 1 )
                         - Tri2D( 2, 1 )
                           ), //====================================================== 2x MEM + OP-a.SUB (REG re-used 2x)
       Tri2D_20_sub_10 = ( Tri2D( 2, 0 )
                         - Tri2D( 1, 0 )
                           ), //====================================================== 2x MEM + OP-b.SUB (REG re-used 2x)
       Tri2D_00_sub_20 = ( Tri2D( 0, 0 )
                         - Tri2D( 2, 0 )
                           ); //====================================================== 2x MEM + OP-c.SUB (REG re-used 1~2x)
//-----------------------
double denominator = ( ( /*
                         Tri2D( 1, 1 )
                       - Tri2D( 2, 1 ) // -------------------------- 2x MEM + OP-1.SUB (avoided by re-use) */
                         Tri2D_11_sub_21 //=========================================== 1x REG + OP-d.MUL
                         ) * ( /*
                               Tri2D( 0, 0 ) //---------------------        + OP-3.MUL
                             - Tri2D( 2, 0 ) //--------------------- 2x MEM + OP-2.SUB (avoided by re-use) */
                               Tri2D_00_sub_20 //===================================== 1x REG + OP-f.ADD
                               ) + ( /*
                                     Tri2D( 2, 0 ) //---------------        + OP-7.ADD
                                   - Tri2D( 1, 0 ) //--------------- 2x MEM + OP-4.SUB (avoided by re-use) */
                                     Tri2D_20_sub_10 //=============================== 1x REG + OP-e.MUL
                                     ) * ( Tri2D( 0, 1 ) //---------        + OP-6.MUL
                                         - Tri2D( 2, 1 ) //--------- 2x MEM + OP-5.SUB
                                           )
                       );
//\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/
//
// Computing the first coefficient ------------------------------------------------------------------------------------------------------
//
double enumer_of_a = ( ( /*
                         Tri2D( 1, 1 )
                       - Tri2D( 2, 1 )  //-------------------------- 2x MEM + OP-8.SUB (avoided by re-use) */
                         Tri2D_11_sub_21 //=========================================== 1x REG + OP-g.MUL
                         ) * ( Point(0)   //------------------------------------------        + OP-i.MUL
                             - Tri2D( 2, 0 ) //--------------------------------------- 2x MEM + OP-h.SUB
                               ) + ( /*
                                     Tri2D( 2, 0 ) //---------------        + OP-E.ADD
                                   - Tri2D( 1, 0 ) //--------------- 2x MEM + OP-B.SUB (avoided by re-use) */
                                     Tri2D_20_sub_10 //=============================== 1x REG + OP-l.ADD
                                     ) * ( Point(1) //--------------------------------        + OP-k.MUL
                                         - Tri2D( 2, 1 ) //--------------------------- 2x MEM + OP-j.MUL
                                           )
                       );/*denominator; *///------------------------ 1x REG + OP-F.DIV (avoided by DEFERRAL THE MOST EXPENSIVE DIVISION UNTIL a third coeff is indeed first needed, if ever-----------[3]

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~J.I.T./non-MISRA-C-RET-->
// TEST CONDITIONS FOR A CHEAPEST EVER J.I.T./non-MISRA-C-RET-->
//
if (  enumer_of_a > denominator                          // in a > 1, THE SIZE DECIDES, the a / denominator > 1, iff enumer_of_a > denominator         a rather expensive .FDIV is avoided at all
   || enumer_of_a * denominator < 0 ) return(false);     // in a < 0, THE SIGN DECIDES, not the VALUE matters, so will use a cheaper .FMUL, instead of a rather expensive .FDIV ~~ ( sign( a ) * sign( denominator ) ) < 0

//\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/
//
// Computing the second coefficient
//
double enumer_of_b = ( ( Tri2D( 2, 1 ) - Tri2D( 0, 1 ) ) //---------------------------------------- 2x MEM + OP-m.SUB
                     * (      Point(0) - Tri2D( 2, 0 ) ) //---------------------------------------- 2x MEM + OP-n.SUB + OP-o.MUL
                     + ( /*
                         Tri2D( 0, 0 ) - Tri2D( 2, 0 )   //--------- 2x MEM + OP-19.SUB + OP-22.ADD (avoided by re-use) */
                         Tri2D_00_sub_20 //======================================================== 1x REG + OP-p.ADD
                         )
                     * (      Point(1) - Tri2D( 2, 1 ) ) //---------------------------------------- 2x MEM + OP-r.SUB + OP-q.MUL
                       );/*denominator; *///------------------------ 1x REG + OP-23.DIV (avoided by DEFERRAL THE MOST EXPENSIVE DIVISION UNTIL a third coeff is indeed first needed, if ever-----------[3]

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~J.I.T./non-MISRA-C-RET-->
// TEST CONDITIONS FOR A 2nd CHEAPEST J.I.T./non-MISRA-C-RET-->
//
if (  enumer_of_b > denominator                          // in b > 1, THE SIZE DECIDES, the a / denominator > 1, iff enumer_of_a > denominator         a rather expensive .FDIV is avoided at all
   || enumer_of_b * denominator < 0 ) return(false);     // in b < 0, THE SIGN DECIDES, not the VALUE matters, so will use a cheaper .FMUL, instead of a rather expensive .FDIV ~~ ( sign( a ) * sign( denominator ) ) < 0

//\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/
//
// Computing the third coefficient
//
double c = 1 - ( ( enumer_of_a
                 - enumer_of_b
                   )
                 / denominator
                 ); // --------------------------------------------- 3x REG + OP-s.SUB + OP-t.FDIC + OP-u.SUB <----THE MOST EXPENSIVE .FDIV BUT HERE, IFF INDEED FIRST NEEDED <---HERE <------------[3]

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~J.I.T./non-MISRA-C-RET-->
// TEST CONDITIONS FOR PRE-FINAL RET J.I.T./non-MISRA-C-RET-->
//
if (   c < 0 || c > 1 ) return( false );

return( true ); //~~~~~~~~~~~~~~~~~~~ "the-last-resort" RET-->

b)算法的其他方法：

让我们回顾一下另一种方法，由于数据更少，指令更少，并且还有望具有高密度，一旦智能使用矢量化 AVX-2 或更好的 AVX-512 向量指令将得到利用每个核心：动画，完全交互式的解释在这里完全与分析问题重新制定。

点到线距离的三重测试（每条线ax + by + c = 0）的廉价成本约为 2FMA3就足够了 + 符号测试（如果 vector-4-in-1-compact AVX-2 / 8- 则更好in-1 AVX-512- VFMADDs)

虽然可能有可能“快速”决定一个点是否有意义，以便针对相应的三角形进行测试，可能polar(R,Theta)通过静态、预先计算的( R_min, R_max, Theta_min, Theta_max )并且“快速”区分每个点，如果它不适合这样的极线段。然而，这样做的成本（数据访问模式成本 + 这些“快速”指令的成本）将超过任何可能“保存”的指令路径，这不需要发生（如果在极线段之外找到一个点）。在 3.0+ GHz 时，每 6 个 CPU 内核约 9 个 CPU 指令的成本实现了 24 个三角形点测试的一系列性能，极地段“预测试”将突然变得异常昂贵，不是说大约二阶负面影响（通过更糟糕的缓存命中/缓存未命中率引入，因为要存储更多数据并将其读入“快速”预测试〜每个三角形“框架”+16B

这显然不是任何进一步行动的好方向，因为性能将下降而不是提高，这是我们在这里的全球战略。

英特尔 i5 8500 CPU可以使用但 AVX-2，因此如果需要，每个内核每个 CPU 时钟节拍使用 8 个三角形的最紧凑的使用方式可以实现 2 倍更高的性能。

TRIPLE-"point-above-line"-TEST per POINT per TRIANGLE:
---------------------------------------------------------------------------------
PRE-COMPUTE STATIC per TRIANGLE CONSTANTS:
                     LINE_1: C0__L1, C1__L1, C2__L1, bool_L1DistanceMustBePOSITIVE
                     LINE_2: C0__L2, C1__L2, C2__L2, bool_L2DistanceMustBePOSITIVE
                     LINE_3: C0__L3, C1__L3, C2__L3, bool_L3DistanceMustBePOSITIVE

TEST per TRIANGLE per POINT (Px,Py) - best executed in an AVX-vectorised fashion
                     LINE_1_______________________________________________________________
                     C0__L1 ( == pre-calc'd CONST = c1 / sqrt( a1^2 + b1^2 ) ) //                       
                Px * C1__L1 ( == pre-calc'd CONST = a1 / sqrt( a1^2 + b1^2 ) ) // OP-1.FMA REG-Px,C1__L1,C0__L1
                Py * C2__L1 ( == pre-calc'd CONST = b1 / sqrt( a1^2 + b1^2 ) ) // OP-2.FMA REG-Py,C2__L1, +
           .GT./.LT. 0                                                         // OP-3.SIG
                     LINE_2_______________________________________________________________
                     C0__L2 ( == pre-calc'd CONST = c2 / sqrt( a2^2 + b2^2 ) ) //                       
                Px * C1__L2 ( == pre-calc'd CONST = a2 / sqrt( a2^2 + b2^2 ) ) // OP-4.FMA REG-Px,C1__L2,C0__L2
                Py * C2__L2 ( == pre-calc'd CONST = b2 / sqrt( a2^2 + b2^2 ) ) // OP-5.FMA REG-Py,C2__L2, +
           .GT./.LT. 0                                                         // OP-6.SIG
                     LINE_3_______________________________________________________________
                     C0__L3 ( == pre-calc'd CONST = c3 / sqrt( a3^2 + b3^2 ) ) //                       
                Px * C1__L3 ( == pre-calc'd CONST = a3 / sqrt( a3^2 + b3^2 ) ) // OP-7.FMA REG-Px,C1__L3,C0__L3
                Py * C2__L3 ( == pre-calc'd CONST = b3 / sqrt( a3^2 + b3^2 ) ) // OP-8.FMA REG-Py,C2__L3, +
           .GT./.LT. 0                                                         // OP-9.SIG

( using AVX-2 intrinsics or inlined assembler will deliver highest performance due to COMPACT 4-in-1 VFMADDs )
                                             ____________________________________________
                                            |  __________________________________________triangle A: C1__L1 
                                            | |  ________________________________________triangle B: C1__L1
                                            | | |  ______________________________________triangle C: C1__L1
                                            | | | |  ____________________________________triandle D: C1__L1
                                            | | | | |
                                            | | | | |      ______________________________
                                            | | | | |     |  ____________________________triangle A: Px
                                            | | | | |     | |  __________________________triangle B: Px
                                            | | | | |     | | |  ________________________triangle C: Px
                                            | | | | |     | | | |  ______________________triandle D: Px
                                            | | | | |     | | | | |
                                            |1|2|3|4|     | | | | |
                                            | | | | |     |1|2|3|4|      ________________
                                            | | | | |     | | | | |     |  ______________triangle A: C0__L1
                                            | | | | |     | | | | |     | |  ____________triangle B: C0__L1
                                            | | | | |     | | | | |     | | |  __________triangle C: C0__L1
                                            | | | | |     | | | | |     | | | |  ________triandle D: C0__L1
                                            | | | | |     | | | | |     | | | | |
                                            |1|2|3|4|     | | | | |     | | | | |
                                            | | | | |     |1|2|3|4|     | | | | |
                                            | | | | |     | | | | |     |1|2|3|4|
  (__m256d)    __builtin_ia32_vfmaddpd256 ( (__v4df )__A, (__v4df )__B, (__v4df )__C ) ~/ per CPU-core @ 3.0 GHz ( for actual uops durations check Agner or Intel CPU documentation )
  can
  perform 4-( point-in-triangle ) PARALLEL-test in just about ~ 9 ASSEMBLY INSTRUCTIONS / per CPU-core @ 3.0 GHz
         24-( point-in-triangle ) PARALLEL-test in just about ~ 9 ASSEMBLY INSTRUCTIONS / per CPU

  using AVX-512 empowered CPU, can use 8-in-1 VFMADDs
  could
  perform 8-( point-in-triangle ) PARALLEL-test in just about ~ 9 ASSEMBLY INSTRUCTIONS / per CPU-core @ 3.0 GHz
         48-( point-in-triangle ) PARALLEL-test in just about ~ 9 ASSEMBLY INSTRUCTIONS / per CPU

c) 最佳并行代码 + 终极性能的提示：

步骤 -1： GPU / CUDA 成本与收益

如果你的博士生导师、教授、老板或项目经理确实坚持要你为这个问题开发一个 GPU 计算的 c++/CUDA 代码解决方案，那么最好的下一步是要求获得任何更适合此类任务的 GPU 卡，而不是您发布的那个。

您指定的卡，Q400 GPU 只有 2 个 SMX（每个 48KB L1 缓存）不太适合真正意义上的并行 CUDA 计算，实际上只有大约 30X 的处理设备来实际执行任何 SIMD 线程块计算，而不是提到它的小内存和微小的 L1 / L2 on-SMX 缓存。因此，在所有与 CUDA 相关的设计和优化工作之后，SIMD 线程中不会有更多，而是一对（！）GPU-SMXwarp32范围内的线程块执行，所以没有大的马戏团可做预计这款基于 GP107 的设备将配备 30+X 更好的设备，用于一些确实高性能的并行处理可用 COTS）

步骤 0：预计算和预排列数据（MAXIMIZE cache-line COHERENCY）：

在这里，优化算法的最佳宏观循环是有意义的，这样您就可以从缓存命中中“受益”最多（即最好地重用“快速”数据，已经预取）

因此，可以测试一下，将工作分配给 N 个并发工作人员是否更快地进行缓存，这些工作人员处理三角形的分离工作块，每个工作块在最小的内存区域上循环——所有点 ( ~ 10k * 2 * 4B ~ 80 kB )，然后进入工作块中的下一个三角形。确保行优先的数组对齐到内存区域是至关重要的（FORTRAN 的家伙可以告诉很多关于这个技巧的成本/收益，对于 HPC 快速紧凑/矢量化向量代数和矩阵处理）

好处？

缓存的系数将被重复使用约 10k 次（成本为~ 0.5~1.0 [ns]，而不是重新获取+ 100 ~ 300 [ns]这些是否必须通过 RAM 内存访问重新读取的成本）。差异，大约~ 200X ~ 600X值得努力使从属于数据访问模式和缓存行资源的工作流最佳地对齐。

结果是原子的 - 任何点都属于一个且仅属于一个（非重叠）三角形。

这简化了对结果向量的先验非冲突和相对稀疏写入，其中任何检测到的三角形内点都可以自由地报告找到的此类三角形的索引。

使用该结果向量可能避免对已经执行匹配（索引为非负）的点进行任何重新测试不是很有效，因为重新读取此类指示和重新安排紧凑对齐的成本4 合 1 或 8 合 1 的三角形点测试，对于不重新测试已经映射的点的潜在“节省”成本将变得不利。

因此，三角形10k块上的点300k可能会产生以下工作流程：
一些300k / 6核心
~ 50k三角形的拆分/每个核心的 1 个核心
~ 50k * 10 k三角形
~ 500M点测试每个核心@ 3.0+ GHz
~ 125M AVX-2向量 4 合 1 紧凑测试 -每个核心的执行次数
~ 125M测试 x10 uops指令@ 3.0 GHz......
这是1.25G uops @ 3.0+ GHz ......秒？是的，有一个非常强烈的动机去这里，朝着这个最终的表现，并以这种方式指导进一步的工作。

所以，我们在这里：

仅在 6 核 AVX-2 i5 8500 @ 3.0+ GHz 上，主要可实现的 HPC 目标在几秒钟内即可实现 300+k 个三角形/10+k 个点

值得努力，不是吗？

于 2019-07-09T16:37:32.700 回答
1

如果您必须执行大量单个查询，四叉树很好，但如果您有 10k 一次要执行所有操作，则有一个专门为此而构建的算法：扫描线。我在单个线程中对类似数据的查询时间不到 5 秒。

我们将画一条垂直线并从左到右扫过您的二维平面，我们将随时跟踪该线与哪些三角形相交。因此，每当您的线在扫描时与您的查询点之一相交时，您只需检查垂直线重叠的三角形。没有其他三角形可以包含该点。

在这种情况下，我有点讨厌“扫描线”这个名字，因为它给人的心理形象是平稳地穿越飞机，但事实并非如此。它只是按从左到右的顺序跳到下一个感兴趣的位置。

根据三角形与查询的比率，您还可以通过它们的 y 坐标将重叠的三角形放在线段树或区间树中，但这实际上对我的数据来说更慢（或者我的实现很糟糕。也可能）。不过绝对值得一试，特别是如果你可以推迟平衡树直到需要它。

设置时间只有 3 种：

  • 获取一个指向三角形的指针向量，按 tri 的最小 x 值排序
  • 获取一个指向三角形的指针向量，按 tri 的最大 x 值排序
  • 获取指向按 x 值排序的点的指针向量

然后我们扫一扫：

  • 将 3 个排序列表中的每一个的“当前位置”初始化为 0
  • 在 3 个排序列表中查找当前位置的最低 x 值。这将是扫描线的下一个 x 位置
    • 如果它是三角形的最左边的点，则将该三角形添加到您的“重叠”向量中
    • 如果它是一个查询点，则针对所有当前重叠的 tris 检查该点
    • 如果它是三角形的最右边点，则将其从重叠向量中删除。
  • 增加该列表的当前位置
  • 循环，直到你的分数，或出三分。

为什么这比四叉树好？因为我们一直在跟踪当前与我们的扫描线相交的 tris。这意味着我们在查询之间重复使用数据，而四叉树并没有真正做到这一点。四叉树将具有更好的单查询性能，但这种批量查找是针对扫描线进行的。

于 2020-02-04T00:56:47.590 回答
1

我认为您可以为每个三角形制作一个带有边界的数组：顶部、底部、右侧和左侧极端。然后，将您的观点与这些界限进行比较。如果它落在一个内，那么看看它是否真的在三角形内。这样，99.9% 的情况不涉及双浮点乘法和一些加法 - 只是比较。仅当该点位于三角形的直线极端内时，才进行计算量大的操作。

这还可以进一步加快，例如通过对三角形的最上端进行排序，并使用二分搜索；然后首先找到三角形下方的最高点，然后检查其上方的点。这样，您只需要检查一半多一点。如果三角形极端的高度有上限，那么您可以检查的次数要少得多。请注意，后一种策略会使您的源代码更加复杂——因此这将是一个确定您希望为优化结果付出多少努力的案例。第一部分似乎相当容易，并且有很大帮助。排序后的列表：更多的努力只是将你的操作减半。我会先看看第一个策略对你来说是否足够。

于 2019-07-08T15:15:44.617 回答
0

使用 boost 的二次 rtree 来找到最近的三角形完美地完成了这项工作。算法现在运行不到一分钟（对于 75k 三角形和 100k 点（10 倍以上！））

我最终通过在每个三角形中放入一个盒子来构建一棵树，并寻找一个点的 Knn 并测试相应的三角形 :) 期待进入像 Spatial 数据库这样的新域会遇到更多麻烦，但 Boost 确实是一个疯狂的库

于 2019-07-09T15:34:22.090 回答