有效查找高密度区域的最佳方法

Posted 2023-02-23

【中文标题】有效查找高密度区域的最佳方法

【英文标题】：Best way to efficiently find high density regions

【发布时间】：2013-11-01 05:43:15

【问题描述】：

在我的编码过程中，我遇到了如下问题： 在二维空间中找到具有最高粒子密度的固定大小的区域。粒子一般可以认为是随机分布在整个空间中，但理论上应该有一些区域密度较高。

例如，100 个粒子随机放置在 500x500 的 2D 网格中，我需要找到 50x50 的区域，其中粒子最多（密度最高）。

除了蛮力测试每个可能的区域（在这种情况下大约超过 200000 个区域）之外，还有其他方法来计算最佳区域吗？对于 n 长度的轴，这将在 O(n^2) 时按比例放大。

【问题讨论】：

@MillieSmith：这是个好主意，但除非您知道需要多少查询，否则它不是很有帮助。

这和这个问题完全一样：poj.org/problem?id=1916 最好的解决方案是O(N^2 + K*D^2) 其中K 是粒子数（在你的情况下是100），D 是区域的大小（在您的情况下为 50），如以下 C++ 代码所示：tausiq.wordpress.com/2012/08/24/uva-10360-rat-attack

【参考方案1】：

算法 1

创建一个 500x500 二维数组，其中每个单元格包含该单元格中粒子数的计数。然后将该数组与 50x50 内核卷积，生成的数组将包含每个单元格中 50x50 区域中的粒子数。然后找到具有最大值的单元格。

如果您使用 50x50 的框作为区域，内核可以分解为两个独立的卷积，每个轴一个。生成的算法是 O(n^2) 空间和时间，其中 n 是您正在搜索的 2D 空间的宽度和高度。

提醒一下，一个带有 boxcar 函数的一维卷积可以在 O(n) 的时间和空间内完成，并且可以就地完成。设 x(t) 为 t=1..n 的输入，设 y(t) 为输出。为 tn 定义 x(t)=0 和 y(t)=0。将内核 f(t) 定义为 0..d-1 的 1 和其他地方的 0。卷积的定义给了我们以下公式：

y(t) = sum i x(t-i) * f(i) = sum i=0..d-1 x(t-i)

这看起来需要时间 O(n*d)，但我们可以将其重写为循环：

y(t) = y(t-1) + x(t) - x(t-d)

这说明一维卷积是O(n)，与d无关。要执行二维卷积，您只需对每个轴执行一维卷积。之所以可行，是因为可以分解 boxcar 内核：通常，大多数内核无法分解。高斯核是另一个可以分解的核，这也是图像编辑程序中的高斯模糊速度如此之快的原因。

对于您指定的数字类型，这将非常快。 500x500 是一个极小的数据集，你的计算机最多可以在几毫秒内检查 202,500 个区域。您将不得不问自己，是否值得额外花费数小时、数天或数周的时间来进一步优化。

这和justhalf的解法一样，只是由于分解了卷积，所以区域大小不影响算法的速度。

算法 2

假设至少有一个点。不失一般性，将二维空间视为整个平面。令 d 为区域的宽度和高度。设 N 为点数。

引理：存在一个密度最大的区域，其左边缘有一个点。

证明：令 R 为最大密度区域。设 R' 是同一个区域，向右平移 R 的左边缘和 R 中最左边的点之间的距离。R 中的所有点也必须位于 R' 中，因此 R' 也是最大密度的区域。

算法

将所有点插入到 K-D 树中。这可以在 O(N log² N) 时间内完成。

对于每个点，考虑宽度为 d 和高度为 2d 的区域，其中该点位于该区域的左边缘。将此区域称为 R。

在 KD 树中查询区域 R 中的点。将此集合称为 S。这可以在 O(N^1/2+|S|) 时间内完成。

找到 R 的 dxd 子区域，其中包含 S 中最大数量的点。这可以在 O(|S| log |S|) 时间内完成，方法是按 y 坐标对 S 进行排序，然后执行线性扫描.

生成的算法的时间为 O(N^3/2 + N |S​​| log |S|)。

比较

当密度高时，算法#1 优于算法#2。算法#2 仅在粒子密度非常低的情况下具有优势，而算法#2 的优势密度会随着总板尺寸的增加而降低。

请注意，可以认为连续情况的密度为零，此时只有算法 #2 有效。

【讨论】：

我不确定你所说的重叠区域是什么意思...我们只搜索一个区域，有 202,500 种可能的检查。

@MillieSmith：性能和优化很重要，但还为时过早。毕竟，优化的第一条规则是“不要”。

@MillieSmith：我提出的方法不依赖于区域的大小。 O() 中的常数因子是非常小的 常数因子。我至少提出了一个完整和正确的解决方案，如果你声称存在一个比无论如何都启发我更有效的解决方案。否则你只是在猜测我错了。

@MillieSmith：我不是在扮演魔鬼的拥护者。我在评论不完整的解决方案。如果它不能解决问题，你的代码有多快也没关系。

@j_random_hacker：分解的框卷积需要每个单元的摊销常数时间。对 d 没有依赖关系，对维数只有线性依赖关系。在一维情况下，给定输入 x(1)..x(n)，输出 y(t) 的公式可以由 y(t) = sum i=0..d-1 x(i +t)，这确实是 O(n*d)。但是，这等价于递归关系 y(t) = y(t-1) + x(t+d-1) - x(t-1)，即 O(n)。这条捷径就是为什么 boxcar 函数的卷积对于离散输入如此有效的原因。

【参考方案2】：

我不知道你使用什么蛮力方法，但最蛮力的方法是O(n^2 d^2)，通过在O(n^2) 时间内迭代每个区域，然后在@987654325 中计算该区域中的粒子数@time 其中d 是您所在区域的大小。

这个问题和这个问题完全一样：Rat Attack，由于区域面积是固定的，所以密度和计数一样， 解决方案是O(n^2 + k*d^2)，其中

n是整个区域的大小（边长） k 是粒子数 d是每个区域的大小（边长）

通过这个算法：

对于每个粒子，更新受此粒子影响的O(d^2) 区域的计数 遍历所有O(n^2)可能的区域，找到最大值

如this code所示，我把相关部分复制到这里供大家参考：

using namespace std;

int mat [1024 + 3] [1024 + 3]; // Here n is assumed to be 1024

int main ()

    int testCases; scanf ("%d", &testCases);

    while ( testCases-- ) 

        Set(mat, 0);

        int d; scanf ("%d", &d); // d is the size of the region
        int k; scanf ("%d", &k); // k is the number of particles

        int x, y, cost;

        for ( int i = 0; i < k; i++ ) 
            scanf ("%d %d %d", &x, &y, &cost); // Read each particle position

            // Update the count of the d^2 region affected by this particle
            for ( int j = max (0, x - d); j <= min (x + d, 1024); j++ ) 
                for ( int k = max (0, y - d); k <= min (y + d, 1024); k++ ) mat [j] [k] += cost;
            
        

        int resX, resY, maxi = -1;

        // Find the maximum count over all regions
        for ( int i = 0; i < 1025; i++ ) 
            for ( int j = 0; j < 1025; j++ ) 
                if ( maxi < mat [i] [j] ) 
                    maxi = mat [i] [j];
                    resX = i;
                    resY = j;
                
            
        

        printf ("%d %d %d\n", resX, resY, maxi);

    
    return 0;

我已将我的 cmets 放入代码中向您解释。

【讨论】：

O(n^2 + k) 的运行时间只需稍作修改即可，因为可以通过 O(n^2) 中的分解框卷积一次完成每个 d*d 大小区域的更新) 时间。

哇，你知道的似乎比我多得多。你能描述一下什么是分解框卷积吗？您还可以在答案中提供一些指向 K-D 树的链接吗？

卷积是获取每个单元格并将其内容添加到相邻单元格的过程。 （当您将值添加到相邻单元格时，您也可以按常数缩放值。）您可以同时对所有单元格执行此操作。内核是您影响的一组相邻单元格，以及每个单元格的系数。 box 函数是一个内核，它只是一个矩形，里面到处都有一个常数值。盒子函数可以分解成两个一维卷积（高斯核也有这个性质），这使得计算卷积非常快。

高斯模糊是卷积的另一个例子，与高斯核的卷积。

【参考方案3】：

将区域划分为 1000x1000 并计算每个（重叠）2x2 中的粒子数。您可以简单地通过规范化 0..1、缩放 0..999 和转换为整数来对它们进行分区。计数可以很容易地存储为整数的二维数组（ushort、uint 或 ulong...mmmm tea）。这相当于在宽相位碰撞检测中使用的简单二维空间分割。