SQRT decomposition 2D

Есть площадь заданного размера, на ней указаны N точек, и даны их координаты. Нужно научиться быстро отвечать на запрос сколько точек находится в квадрате. Моя идея заключается в том, чтобы поделить область на sqrt(N) * sqrt(N) меньших областей, преподсчетом решить задачу на них, по данной матрице построить матрицу частичных сумм, и в итоге просто после поступления запроса (координат квадрата) добавлять к подсчитанным точкам из маленьких квадратов, рассмотренные отдельно, вручную. Что-то такое:

int sqrt_count = std::ceil(std::sqrt(points_num));
int rect_size = std::ceil(lim / sqrt_count);
matrix2d ppr = matrix2d(sqrt_count, std::vector<int>(sqrt_count, 0)); // points_per_rect
for (int i = 0; i < points_num; ++i)
{
    ++ppr[points[i].X / rect_size][points[i].Y / rect_size];
}
matrix2d blocks = matrix2d(sqrt_count + 1, std::vector<int>(sqrt_count + 1, 0));
for (int i = 1; i < sqrt_count; ++i)
{
    for (int j = 1; j < sqrt_count; ++j)
    {
        blocks[i][j] = blocks[i][j - 1] + blocks[i - 1][j] - blocks[i - 1][j - 1] + ppr[i - 1][j - 1];
    }
}

Но проблема с которой я столкнулся, это то что я не пойму, как теперь мне узнать, какие конкретно точки принадлежат какому из этих квадратов. Допустим поступает запрос, где используются два цельных маленьких квадрата, и один наполовину. Два цельных я могу просто добавить не считая, а для того что наполовину, мне нужно проверить какие точки из него входят, а какие нет. Хранить отдельно те точки что я добавляю для каждого квадратика очень накладно. В обычной одномерной sqrt декомпозиции все просто, ведь там преобразовывается индекс, а как быть здесь?