更好的计算网格点数的算法

better algorithm for counting points in grids

本文关键字：算法 网格 计算 更好      更新时间：2023-10-16

假设我在[0,1]*[0,1]中有一个点p，并且[0,1]被划分为m(假设为200)个网格。我使用A[m][m]来表示[以P为中心，长度为2h的小正方形]是否覆盖每个网格。因此对于点P, A[i][j]为(递增)1或0。

假设我有n个这样的点(P1，…，Pn)，我想计算A(对于每个点Pi，我重做上述过程，加1或不加1)。我如何有效地做到这一点(用c++)，而不是编写3层for循环来检查每个网格和每个点(所以O(nm^2))?

我用c++尝试了简单的3 for循环。它比在r中使用一些矢量化操作(如vector<= number用于比较n个数，A[bool vector, bool vector]用于子集)需要更长的时间。

由于c++通常比R更快，有什么聪明的方法来实现这个过程吗?

#include <Rcpp.h>
#include <cmath>
using namespace Rcpp;
// [[Rcpp::export]]
double myfun(NumericVector u, NumericVector v)
{
    double n = u.size();
    double A[200][200] = {0};
    double pos[200];
    int i = 0, j = 0, k = 0;
    for (i = 0; i < 200; i++)
    {
        pos[i] = (double)i / 201;
    }
    for (k = 0; k < n; k++)
    {
        for (i = 0; i < 200; i++)
        {
            for (j = 0; j < 200; j++)
            {
                if ( (fabs(u[k] - pos[i]) <= h) && (fabs(v[k] - pos[j]) <=h ) )
                {
                    A[i][j]++;
                }
            }
        }
    }
    double s = 0, avg = 0;
    for (i = 0; i <200; i++)
    {
        for (j = 0; j < 200; j++)
        {
            s += A[i][j];
        }
    }
    avg = s / (200 * 200);
    return (avg);
}

这两个内部循环只决定网格中点的索引。但是你可以直接计算索引:

int i = (int)(u[k]*200);
int j = (int)(v[k]*200);

你可能还需要检查i和j没有达到索引200。这只发生在u[k] == 1.0或v[k] == 1.0。

double n = u.size();
double A[200][200] = {0};
for (int k = 0; k < n; k++)
{
    int i = (int)(u[k]*200);
    int j = (int)(v[k]*200);
    if (i == 200)
        i = 199;
    if (j == 200)
        j = 199;
    A[i][j]++;
}
相关文章：