高效算法获取圆心周围的点

Question

高效算法获取圆心周围的点

c++algorithmcanvasgeometrydistance

10

问题

我想要获取给定半径内，以给定点为圆心的所有像素，其中点只能具有整数坐标，即画布上的像素。

我想获取给定 (x, y) 和 r 的黄色区域内的所有点。

方法

我能想到的最有效的方法是循环遍历围绕 (x, y) 的正方形，并检查每个点的 欧几里得距离：

for (int px = x - r; px <= x + r; px++) {
  for (int py = y - r; py <= y + r; py++) {
    int dx = x - px, dy = y - py;

    if (dx * dx + dy * dy <= r * r) {
      // Point is part of the circle.
    }
  }
}

然而，这意味着该算法将检查不属于圆形的(r * 2)^2 * (4 - pi) / 4像素。dx * dx + dy * dy <= r * r似乎相当昂贵，被冗余地称为几乎1/4的时间。

集成类似此处提出的内容可能会提高性能。

for (int px = x - r; px <= x + r; px++) {
  for (int py = y - r; py <= y + r; py++) {
    int dx = abs(x - px), dy = abs(y - py);

    if (dx + dy <= r || (!(dx > r || dy > r) && (dx * dx + dy * dy <= r * r))) {
      // Point is part of the circle.
    }
  }
}

然而，正如作者所指出的那样，在大多数点都在圆内时（特别是由于 abs），这很可能不会更快，而在这种情况下，pi / 4 就在其中。

我没有找到关于这个问题的任何资源。我特别寻求C++解决方案，而不是SQL相关的内容。

- creativecreatorormaybenot

1

整数乘法和加法非常快。我不确定添加大量额外的比较会有所帮助，特别是因为它可能引入所有可能的分支预测错误。 - François Andrieux

通常情况下，“最好的”或“最佳的”问题往往无法得到明确的答案，因为很少有单一的最佳解决方案。哪种解决方案最好取决于实际使用情况和它将运行的架构。 - François Andrieux

3

利用已知给定y值的x值范围可以帮助你了解前一个和后一个y值对应的x值范围，这是一个优势。 - Scott Hunter

1

一旦你从前一个点变坏到下一个点，你至少应该能够从循环中跳出。此时，你就知道已经越过了边界，在该线上没有更多的好点了。 - NathanOliver

1

另一个选项是，为了制作Scotts的示例，从中间的顶部或底部开始，然后循环到左边和右边，直到你碰到边界。这意味着当你通过边界时，内部循环会立即停止。 - NathanOliver

我认为二分查找可以改进解决方案。 - Mandy007

7个回答

3

代码

基于@ScottHunter的想法，我提出了以下算法：

#include <functional>

// Executes point_callback for every point that is part of the circle
// defined by the center (x, y) and radius r.
void walk_circle(int x, int y, int r,
                 std::function<void(int x, int y)> point_callback) {
  for (int px = x - r; px < x + r; px++)
    point_callback(px, y);
  int mdx = r;
  for (int dy = 1; dy <= r; dy++)
    for (int dx = mdx; dx >= 0; dx--) {
      if (dx * dx + dy * dy > r * r)
        continue;
      for (int px = x - dx; px <= x + dx; px++) {
        point_callback(px, y + dy);
        point_callback(px, y - dy);
      }
      mdx = dx;
      break;
    }
}

算法解释

该算法执行了少量的检查。具体来说，它仅在每行中检查直到第一个属于圆的点被找到。此外，在下一行中，它将跳过先前已识别点左侧的点。此外，通过使用对称性，只检查了一半的行（从0开始为n/2 + 1/2）。

可视化

这是我创建的算法的可视化。红色轮廓表示以前将要检查的正方形，黑色像素表示实际圆（红色像素在中心）。算法检查点（标记为蓝色），并循环遍历有效点（标记为绿色）。
正如您所看到的，最后的蓝色像素数量是微不足道的，即只有很少的点被循环遍历而不是圆的一部分。此外，请注意，每次只需要检查第一个绿色像素，其他像素只需循环遍历，这就是为什么它们会立即出现的原因。

注释

轴可以很容易地反转，显然。

这可以通过更多地利用对称性进行优化，即行将与列相同（通过所有行与从左到右，从上到下，反之亦然的所有列相同），从中心向下仅遍历四分之一的行即可确定哪些点将成为圆的一部分。但是，我觉得这样做带来的轻微性能提升不值得额外的代码。
如果有人想编写它，请提出对此答案的修改建议。

带注释的代码

#include <functional>

// Executes point_callback for every point that is part of the circle
// defined by the center (x, y) and radius r.
void walk_circle(int x, int y, int r,
                 std::function<void(int x, int y)> point_callback) {
  // Walk through the whole center line as it will always be completely
  // part of the circle.
  for (int px = x - r; px < x + r; px++)
    point_callback(px, y);
  // Define a maximum delta x that shrinks whith every row as the arc
  // is closing.
  int mdx = r;
  // Start directly below the center row to make use of symmetry.
  for (int dy = 1; dy <= r; dy++)
    for (int dx = mdx; dx >= 0; dx--) {
      // Check if the point is part of the circle using Euclidean distance.
      if (dx * dx + dy * dy > r * r)
        continue;

      // If a point in a row left to the center is part of the circle,
      // all points to the right of it until the center are going to be
      // part of the circle as well.
      // Then, we can use horizontal symmetry to move the same distance
      // to the right from the center.
      for (int px = x - dx; px <= x + dx; px++) {
        // Use y - dy and y + dy thanks to vertical symmetry
        point_callback(px, y + dy);
        point_callback(px, y - dy);
      }

      // The next row will never have a point in the circle further left.
      mdx = dx;
      break;
    }
}

- creativecreatorormaybenot

JS 可能足够快，可以“暴力计算”每个圆的像素坐标（可能不需要优化）。但是，您只测试一个半径或仅测试几个半径吗？如果是，您可以将圆的内部点预先计算到数组中并重复使用这些计算出的点 - 偏移到您正在测试的圆的原点。 - markE

1

你可以看到我的代码的最后更新，我认为那是更好的。 - Mandy007

2

该问题的固定复杂度为O（n^2），其中n是圆的半径。与正方形或任何常规的2D形状相同的复杂度。

无法否认的事实是，即使利用对称性，也不能减少圆中像素的数量，因此复杂度仍然相同。

因此，忽略复杂性并寻求优化。

在您的问题中，您声明每个像素（或第4个像素）中的abs有点太昂贵。

每行一次比每个像素一次更好

您可以将其减少到每行1个平方根。对于半径256的圆来说，这是128个平方根。

void circle(int x, int y, int radius) {
    int y1 = y, y2 = y + 1, r = 0, rSqr = radius * radius;
    while (r < radius) {
        int x1 = x, x2 = x + 1, right = x + sqrt(rSqr - r * r) + 1.5;;
        while (x2 < right) {
            pixel(x1, y1);
            pixel(x2, y1);
            pixel(x1--, y2);
            pixel(x2++, y2);
        }
        y1--;
        y2++;
        r++;
    }
}

为了更好地利用它，您可以创建一个查找表来进行sqrt根计算。

所有整数

或者，您可以使用变体的Bresenham线算法，用所有整数数学代替平方根。但是这很混乱，并且除非设备没有浮点单元，否则不会有任何好处。

void circle(int x, int y, int radius) {
    int l, yy = 0, xx = radius - 1, dx = 1, dy = 1;
    int err = dx - (radius << 1);
    int l2 = x, y0 = y, r2 = x + 1;
    int l1 = x - xx, r1 = r2 + xx;
    int y2 = y0 - xx, y1 = y0 + 1, y3 = y1 + xx;
    while (xx >= yy) {
        l = l1;
        while (l < r1) {
            pixel(l, y1);
            pixel(l++, y0);
        }
        l = l2;
        while (l < r2) {
            pixel(l, y3);
            pixel(l++, y2);
        }
        err += dy;
        dy += 2;
        y0--;
        yy++;
        y1++;
        l2--;
        r2++;
        if (err > 0) {
            dx += 2;
            err += (-radius << 1) + dx;
            xx--;
            r1--
            l1++
            y3--
            y2++
        }
    }
}

- Blindman67

2

这是一种优化方法，可以将搜索的维度减少1/4：

for (int px = x; px <= x + r; ++px) {
  bool find = false;
  int dx = x - px, dy;
  for (int py = y; !find && py <= y + r; ++py) {
    dy = y - py;
    if (dx * dx + dy * dy <= r * r)) {
      /* (px, py), (px, y+y-py+r), (x+x-px+r, py) 
       & (x+x-px+r, y+y-py+r) are part of the circle.*/
    }else{
      find = true; //Avoid increasing on the axis y
    }
  }
}

更好的方式是通过优化第二个循环中的迭代，避免使用if条件语句来提高性能。

for (int px = x; px <= x + r; ++px) {
  int dx = x - px, py = y;
  for (; dx * dx + (py-y) * (py-y) <= r * r; ++py) {
    /* (px, py), (px, y+y-py+r), (x+x-px+r, py) 
     & (x+x-px+r, y+y-py+r) are part of the circle.*/
  }
}

我认为另一个选项是使用二分查找来寻找上限：

int binarySearch(int R, int dx, int y){
  int l=y, r=y+R;
  while (l < r) { 
    int m = l + (r - l) / 2;  
    if(dx*dx + (y - m)*(y - m) > R*R) r = m - 1; 
    else if(dx*dx + (y - m)*(y - m) < R*R) l = m + 1; 
    else r = m;
  }
  if(dx*dx + (y - l)*(y - l) > R*R) --l;
  return l;
}

for (int px = x; px <= x + r; ++px) {
  int upperLimit = binarySearch(r, px-x, y);
  for (int py = y; py <= upperLimit; ++py) {
    /* (px, py), (px, y+y-py+r), (x+x-px+r, py) 
     & (x+x-px+r, y+y-py+r) are part of the circle.*/
  }
}

二分查找的思想是要最优地找到上限，避免在 for 循环中进行 if 条件判断和计算。为此，需要检查哪个整数最大，使得当前点与圆周半径之间的距离最小。

注：抱歉我的英语不好。

- Mandy007

2

for (line = 1; line <= r; line++) {
   dx = (int) sqrt(r * r - line * line);
   for (ix = 1; ix <= dx; ix++) {
       putpixel(x - ix, y + line)
       putpixel(x + ix, y + line)
       putpixel(x - ix, y - line)
       putpixel(x + ix, y - line)
   } 
}

为了避免在坐标轴处重复生成像素，值得从1开始循环，并在单独的循环中绘制中心线（ix==0或line==0）。

请注意，还有一种纯整数Bresenham算法用于生成圆周点。

- MBo

这里使用对称性没有真正的好处。我改变了代码以生成一个四分之一圆形。 - MBo

2

首先我们计算圆的内切正方形。这个公式很简单：

x² + y² = r²    // circle formula
2h² = r²        // all sides of square are of equal length so x == y, lets define h := x
h = r / sqrt(2) // half side length of the inner square

现在，圆内的每个点都位于(-h, -h)和(+h, +h)之间。这是我的意思：

剩下的蓝色部分有点棘手，但也不太复杂。我们从蓝色圆的顶部开始(x=0，y=-radius)。接下来，我们向右走(x++)直到离开圆周为止(直到x²+y²

(-x，-y)，(+x，-y) (-x，+y)，(+x，+y) (-y，-x)，(-y，+x) (+y，-x)，(+y，+x)

现在我们向下移动1行(y--)并重复上面的步骤(同时保持最近的x值)。将圆心添加到每个点中，就完成了。

这是一个可视化效果。由于放大缩小的原因，存在一些伪像。红色点显示了我们在每次迭代时测试的内容：

这是完整的代码(使用opencv绘制图形)：

#include <opencv2/opencv.hpp>

constexpr double sqrt2 = 1.41421356237309504880168;

int main()
{
    cv::Point center(200, 200);
    constexpr int radius = 180;

    // create test image
    cv::Mat img(400, 400, CV_8UC3);
    cv::circle(img, center, radius, {180, 0, 0}, cv::FILLED);
    cv::imshow("img", img);
    cv::waitKey();

    // calculate inner rectangle
    int halfSideLen = radius / sqrt2;
    cv::Rect innerRect(center.x - halfSideLen, center.y - halfSideLen, halfSideLen * 2, halfSideLen * 2);
    cv::rectangle(img, innerRect, {0, 180, 0}, cv::FILLED);
    cv::imshow("img", img);
    cv::waitKey();

    // probe the rest
    int x = 0;
    for (int y = radius; y >= halfSideLen; y--)
    {
        for (; x * x + y * y < radius * radius; x++)
        {
            // anything between the following points lies within the circle
            // each pair of points represents a line
            // (-x, -y), (+x, -y)
            // (-x, +y), (+x, +y)
            // (-y, -x), (-y, +x)
            // (+y, -x), (+y, +x)

            // center + {(-X..X) x (-Y..Y)} is inside the circle
            cv::line(img, cv::Point(center.x - x, center.y - y), cv::Point(center.x + x, center.y - y), {180, 180, 0});
            cv::line(img, cv::Point(center.x - x, center.y + y), cv::Point(center.x + x, center.y + y), {180, 180, 0});
            cv::line(img, cv::Point(center.x - y, center.y - x), cv::Point(center.x - y, center.y + x), {180, 180, 0});
            cv::line(img, cv::Point(center.x + y, center.y - x), cv::Point(center.x + y, center.y + x), {180, 180, 0});

            cv::imshow("img", img);
            cv::waitKey(20);
        }
    }

    cv::waitKey();
    return 0;
}

- Timo

@creativecreatorormaybenot 我认为我的版本应该比较快，除了非常小的图像之外，因为中心正方形的计算是恒定时间（没有任何三角函数，我们只使用基本数字）。这也意味着您不必在正方形内执行任何检查。对于图像的其余部分，算法基本相同。 - Timo

1

刚刚也注意到了。这可能会归结为汇编级别的微小优化。也许我稍后会进行基准测试。 - Timo

在运行算法之前，是否有关于像素数量的解析解？ - D Duck

@DDuck 你可以通过计算圆的面积来近似计算像素数。但是，该结果的准确性取决于用于绘制圆的算法。 - Timo

1

你可以画一个适合圆内的正方形，很容易找到点是否在其中。这将在O(1)时间内解决大多数点（2 * r^2），而不是搜索所有点（4 * r^2）。

编辑：对于剩下的点，您不需要遍历所有其他像素。您需要遍历4个矩形，大小为[(2r/sqrt(2)), r-(r/sqrt(2))]，位于内部正方形的4个边缘（北、东、南、西）。这意味着您永远不必搜索角落上的正方形。由于它完全对称，我们可以取输入点的绝对值，并搜索点是否在坐标平面正半轴上的半个正方形中。这意味着我们只需要循环一次，而不是4次。

int square_range = r/sqrt(2);
int abs_x = abs(x);
int abs_y = abs(y);

if(abs_x < square_range && abs_y < square_range){
    //point is in
}
else if(abs_x < r && abs_y < r){  // if it falls in the outer square
    // this is the only loop that has to be done
    if(abs_x < abs_y){
        int temp = abs_y;
        abs_y = abs_x;
        abs_x = temp;
    }
    for(int x = r/sqrt(2) ; x < r ; x++){
        for(int y = 0 ; y < r/sqrt(2) ; y++){
             if(x*x + y*y < r*r){
                 //point is in
             }
         }    
    }    
}

代码的总复杂度为O((r-r/sqrt(2))* (r/sqrt(2)))。这仅循环了内部正方形和圆形外边界之间的一个矩形（8向对称）。

- heapoverflow

这种方法会在圆内画一个正方形。那么剩下的圆呢？ - MBo

正如我之前提到的，超出方框范围的其他点将通过循环进行搜索。 - heapoverflow

整体复杂度是O(n^2)，而不是您提供的O((r-r/sqrt(2))* (r/sqrt(2)))的复杂度。二次方和大O符号不应该包括系数，只有最高幂项才需要考虑。这个问题无法简化到低于O(n^2)的复杂度。 - Blindman67

我知道大O符号表示法，但我的观点是指出这个算法与其他O(n^2)相比执行的操作更少。 - heapoverflow

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Timo · Accepted Answer

好的，这里是我承诺的基准测试结果。

设置

我使用了google benchmark，任务是将圆周上所有点插入到std::vector<point>中。我对一组半径和一个固定的圆心进行基准测试：

radii = {10, 20, 50, 100, 200, 500, 1000}
center = {100, 500}

语言: C++17
编译器: msvc 19.24.28316 x64
平台: Windows 10
优化: O2 (完全优化)
线程处理: 单线程执行

每个算法的结果都会进行正确性测试（与 OP 算法的输出进行比较）。

目前已经对以下算法进行了基准测试：

OP 的算法 enclosing_square。
我的算法 containing_square。
creativecreatorormaybenot 的算法 edge_walking。
Mandy007 的算法 binary_search。

结果

Run on (12 X 3400 MHz CPU s)
CPU Caches:
  L1 Data 32K (x6)
  L1 Instruction 32K (x6)
  L2 Unified 262K (x6)
  L3 Unified 15728K (x1)
-----------------------------------------------------------------------------
Benchmark                                   Time             CPU   Iterations
-----------------------------------------------------------------------------
binary_search/10/manual_time              804 ns         3692 ns       888722
binary_search/20/manual_time             2794 ns        16665 ns       229705
binary_search/50/manual_time            16562 ns       105676 ns        42583
binary_search/100/manual_time           66130 ns       478029 ns        10525
binary_search/200/manual_time          389964 ns      2261971 ns         1796
binary_search/500/manual_time         2286526 ns     15573432 ns          303
binary_search/1000/manual_time        9141874 ns     68384740 ns           77
edge_walking/10/manual_time               703 ns         5492 ns       998536
edge_walking/20/manual_time              2571 ns        49807 ns       263515
edge_walking/50/manual_time             15533 ns       408855 ns        45019
edge_walking/100/manual_time            64500 ns      1794889 ns        10899
edge_walking/200/manual_time           389960 ns      7970151 ns         1784
edge_walking/500/manual_time          2286964 ns     55194805 ns          308
edge_walking/1000/manual_time         9009054 ns    234575321 ns           78
containing_square/10/manual_time          629 ns         4942 ns      1109820
containing_square/20/manual_time         2485 ns        40827 ns       282058
containing_square/50/manual_time        15089 ns       361010 ns        46311
containing_square/100/manual_time       62825 ns      1565343 ns        10990
containing_square/200/manual_time      381614 ns      6788676 ns         1839
containing_square/500/manual_time     2276318 ns     45973558 ns          312
containing_square/1000/manual_time    8886649 ns    196004747 ns           79
enclosing_square/10/manual_time          1056 ns         4045 ns       660499
enclosing_square/20/manual_time          3389 ns        17307 ns       206739
enclosing_square/50/manual_time         18861 ns       106184 ns        37082
enclosing_square/100/manual_time        76254 ns       483317 ns         9246
enclosing_square/200/manual_time       421856 ns      2295571 ns         1654
enclosing_square/500/manual_time      2474404 ns     15625000 ns          284
enclosing_square/1000/manual_time     9728718 ns     68576389 ns           72

代码

下面是完整的测试代码，您可以复制并粘贴它进行测试。 fill_circle.cpp包含了不同算法的实现。

main.cpp

#include <string>
#include <unordered_map>
#include <chrono>

#include <benchmark/benchmark.h>

#include "fill_circle.hpp"

using namespace std::string_literals;

std::unordered_map<const char*, circle_fill_func> bench_tests =
{
    {"enclosing_square", enclosing_square},
    {"containing_square", containing_square},
    {"edge_walking", edge_walking},
    {"binary_search", binary_search},
};

std::vector<int> bench_radii = {10, 20, 50, 100, 200, 500, 1000};

void postprocess(std::vector<point>& points)
{
    std::sort(points.begin(), points.end());
    //points.erase(std::unique(points.begin(), points.end()), points.end());
}

std::vector<point> prepare(int radius)
{
    std::vector<point> vec;
    vec.reserve(10ull * radius * radius);
    return vec;
}

void bm_run(benchmark::State& state, circle_fill_func target, int radius)
{
    using namespace std::chrono;
    constexpr point center = {100, 500};

    auto expected_points = prepare(radius);
    enclosing_square(center, radius, expected_points);
    postprocess(expected_points);

    for (auto _ : state)
    {
        auto points = prepare(radius);

        auto start = high_resolution_clock::now();
        target(center, radius, points);
        auto stop = high_resolution_clock::now();

        postprocess(points);
        if (expected_points != points)
        {
            auto text = "Computation result incorrect. Expected size: " + std::to_string(expected_points.size()) + ". Actual size: " + std::to_string(points.size()) + ".";
            state.SkipWithError(text.c_str());
            break;
        }

        state.SetIterationTime(duration<double>(stop - start).count());
    }
}

int main(int argc, char** argv)
{
    for (auto [name, target] : bench_tests)
        for (int radius : bench_radii)
            benchmark::RegisterBenchmark(name, bm_run, target, radius)->Arg(radius)->UseManualTime();

    benchmark::Initialize(&argc, argv);
    if (benchmark::ReportUnrecognizedArguments(argc, argv))
        return 1;
    benchmark::RunSpecifiedBenchmarks();
}

fill_circle.hpp

#pragma once

#include <vector>

struct point
{
    int x = 0;
    int y = 0;
};

constexpr bool operator<(point const& lhs, point const& rhs) noexcept
{
    return lhs.x != rhs.x
               ? lhs.x < rhs.x
               : lhs.y < rhs.y;
}

constexpr bool operator==(point const& lhs, point const& rhs) noexcept
{
    return lhs.x == rhs.x && lhs.y == rhs.y;
}

using circle_fill_func = void(*)(point const& center, int radius, std::vector<point>& points);

void enclosing_square(point const& center, int radius, std::vector<point>& points);
void containing_square(point const& center, int radius, std::vector<point>& points);
void edge_walking(point const& center, int radius, std::vector<point>& points);
void binary_search(point const& center, int radius, std::vector<point>& points);

fill_circle.cpp

#include "fill_circle.hpp"

constexpr double sqrt2 = 1.41421356237309504880168;
constexpr double pi = 3.141592653589793238462643;

void enclosing_square(point const& center, int radius, std::vector<point>& points)
{
    int sqr_rad = radius * radius;

    for (int px = center.x - radius; px <= center.x + radius; px++)
    {
        for (int py = center.y - radius; py <= center.y + radius; py++)
        {
            int dx = center.x - px, dy = center.y - py;
            if (dx * dx + dy * dy <= sqr_rad)
                points.push_back({px, py});
        }
    }
}

void containing_square(point const& center, int radius, std::vector<point>& points)
{
    int sqr_rad = radius * radius;
    int half_side_len = radius / sqrt2;
    int sq_x_end = center.x + half_side_len;
    int sq_y_end = center.y + half_side_len;

    // handle inner square
    for (int x = center.x - half_side_len; x <= sq_x_end; x++)
        for (int y = center.y - half_side_len; y <= sq_y_end; y++)
            points.push_back({x, y});

    // probe the rest
    int x = 0;
    for (int y = radius; y > half_side_len; y--)
    {
        int x_line1 = center.x - y;
        int x_line2 = center.x + y;
        int y_line1 = center.y - y;
        int y_line2 = center.y + y;

        while (x * x + y * y <= sqr_rad)
            x++;

        for (int i = 1 - x; i < x; i++)
        {
            points.push_back({x_line1, center.y + i});
            points.push_back({x_line2, center.y + i});
            points.push_back({center.x + i, y_line1});
            points.push_back({center.x + i, y_line2});
        }
    }
}

void edge_walking(point const& center, int radius, std::vector<point>& points)
{
    int sqr_rad = radius * radius;
    int mdx = radius;

    for (int dy = 0; dy <= radius; dy++)
    {
        for (int dx = mdx; dx >= 0; dx--)
        {
            if (dx * dx + dy * dy > sqr_rad)
                continue;

            for (int px = center.x - dx; px <= center.x + dx; px++)
            {
                for (int py = center.y - dy; py <= center.y + dy; py += 2 * dy)
                {
                    points.push_back({px, py});
                    if (dy == 0)
                        break;
                }
            }

            mdx = dx;
            break;
        }
    }
}

void binary_search(point const& center, int radius, std::vector<point>& points)
{
    constexpr auto search = []( const int &radius, const int &squad_radius, int dx, const int &y)
    {
        int l = y, r = y + radius, distance;

        while (l < r)
        {
            int m = l + (r - l) / 2;
            distance = dx * dx + (y - m) * (y - m);
            if (distance > squad_radius)
                r = m - 1;
            else if (distance < squad_radius)
                l = m + 1;
            else
                r = m;
        }

        if (dx * dx + (y - l) * (y - l) > squad_radius)
            --l;

        return l;
    };

    int squad_radius = radius * radius;    
    for (int px = center.x - radius; px <= center.x + radius; ++px)
    {
        int upper_limit = search(radius, squad_radius, px - center.x, center.y);
        for (int py = 2*center.y - upper_limit; py <= upper_limit; ++py)
        {
            points.push_back({px, py});
        }
    }
}