Эффективный алгоритм для получения точек в круге вокруг центра - PullRequest
8 голосов
/ 05 марта 2020

Задача

Я хочу получить все пиксели, которые находятся в пределах круга данного радиуса относительно данной точки, где точки могут иметь только целочисленные координаты , то есть пиксели на холсте.

Illustration

Поэтому я хочу получить все точки в желтой области с учетом (x, y) и r.

Подходы

Самый эффективный способ, который я могу придумать, - это l oop через квадрат вокруг (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.

Ответы [ 7 ]

4 голосов
/ 06 марта 2020

Хорошо, вот тесты, которые я обещал.

Настройка

Я использовал google benchmark , и задача состояла в том, чтобы вставить все точки внутри периметра круга в std::vector<point>. Я бенчмарк для набора радиусов и постоянного центра:

radii = {10, 20, 50, 100, 200, 500, 1000}
center = {100, 500}
  • язык: C ++ 17
  • компилятор: MSV c 19.24.28316 x64
  • платформа: windows 10
  • оптимизация: O2 (полная оптимизация)
  • многопоточность: однопоточное выполнение

Результаты каждого алгоритма проверяются на корректность (по сравнению с выводом алгоритма OP).

До сих пор тестировались следующие алгоритмы:

  1. Алгоритм OP enclosing_square.
  2. Мой алгоритм containing_square.
  3. алгоритм creativecreatormamabenot edge_walking.
  4. алгоритм 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});
        }
    }
}
2 голосов
/ 07 марта 2020

Задача имеет фиксированную сложность O (n ^ 2), где n - радиус круга. Сложность та же, что и у квадрата или любой другой двумерной фигуры

Нельзя обойти стороной тот факт, что вы не можете уменьшить количество пикселей в круге, даже если вы используете симметрию, сложность остается прежней. .

Так что игнорируя сложность и ищите оптимизацию.

В своем вопросе вы утверждаете, что abs немного дороже на пиксель (или 4-й пиксель)

Один раз на строку лучше, чем один раз на пиксель

Вы можете уменьшить его до 1 квадрата root на строку. Для радиуса окружности 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 root.

Все целые числа

В качестве альтернативы вы можете использовать вариацию линии Брезенхама , которая заменяет квадрат root всей целой математикой. Однако это беспорядок и не принесет никакой пользы, если устройство не имеет блока с плавающей запятой.

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++
        }
    }
}
2 голосов
/ 05 марта 2020

Хорошо, прежде всего мы вычисляем внутренний квадрат круга. Формула для него проста:

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) лежит внутри круга. Вот изображение того, что я имею в виду:

1

Оставшаяся синяя часть немного хитрая, но и не слишком сложная. Мы начинаем с самой вершины синего круга (x = 0, y = -radius). Затем мы идем направо (x++) до тех пор, пока не покинем периметр круга (пока x²+y² < r² больше не будет удерживаться). Все между (0, y) и (x, y) находится внутри круга. Из-за симметрии мы можем расширить эту 8-кратную величину на

  • (- x, -y), (+ x, -y)
  • (- x, + y), (+ x, + y)
  • (- y, -x), (-y, + x)
  • (+ y, -x), (+ y, + x)

теперь мы go опустились на 1 строку (y--) и повторили шаги, описанные выше (сохраняя самое последнее значение x). Добавьте центр круга к каждой из точек, и все готово.

Вот визуализация. Есть некоторые артефакты из-за апскейлинга. Красная точка показывает то, что мы тестируем на каждой итерации:

1

Вот полный код (используя 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;
}
2 голосов
/ 05 марта 2020
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) в отдельных l oop.

Обратите внимание, что существует также чистый целочисленный алгоритм Брезенхэма для генерации точек окружности.

2 голосов
/ 05 марта 2020

Код

Основываясь на идее @ 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;
    }
}

Алгоритм объяснен

Это Алгоритм выполняет минуту количество проверок. В частности, он проверяет только в каждой строке, пока не будет достигнута первая точка, которая является частью круга. Кроме того, он будет пропускать точки слева от ранее определенной точки в следующем ряду. Кроме того, используя симметрию, проверяется только половина строк (n/2 + 1/2, начиная с 0).

Visualization

Это визуализация созданного мною алгоритма. Красный контур обозначает квадрат, который ранее был проверен, а черные пиксели - реальный круг (красный пиксель в центре является центром). Алгоритм проверяет точки (отмечены синим цветом) и проходит через действительные точки (отмечены зеленым цветом).
Как видите, количество синих пикселей в конце - это минуты, т. Е. Есть только несколько точек, которые зациклены на часть круга. Кроме того, обратите внимание, что только первый зеленый пиксель каждый раз нуждается в проверке, остальные только зацикливаются, поэтому они появляются мгновенно.

Примечания

Оси можно легко поменять местами, очевидно .

Это можно оптимизировать, если использовать симметрию еще больше, т. Е. Строки будут такими же, как столбцы (прохождение всех строк аналогично прохождению всех столбцов слева направо, вверх, вниз, наоборот (vise vera) и спуск только на четверть рядов от центра будет достаточно, чтобы точно определить, какие точки будут частью круга. Однако я чувствую, что незначительное снижение производительности, которое это собирается дать, не стоит дополнительного кода.
Если кто-то хочет его закодировать, предложите изменить этот ответ.

Код с комментариями

#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;
    }
}
2 голосов
/ 05 марта 2020

Это оптимизация, которая уменьшает 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
    }
  }
}

или лучше, улучшая производительность итерации второго круга for, избегая условного 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.*/
  }
}

Идея бинарного поиска состоит в том, чтобы найти верхний предел оптимально, избегая условия if и вычислений в пределах for цикл. Для этого проверяется, какое наибольшее целое число определяет расстояние между текущей точкой и радиусом внутри круга.

PD: Извините, мой английский sh.

1 голос
/ 05 марта 2020

Вы можете нарисовать квадрат, который помещается внутри круга, и довольно просто определить, попадет ли точка в точку.

Это решит большинство точек (2 * r ^ 2) за O (1) время вместо поиска всех (4 * r ^ 2) точек.

Edit: Для остальных точек вам не нужно l oop все остальные пиксели. Вам нужно l oop для 4 прямоугольников с размерами [(2r / sqrt (2)), r- (r / sqrt (2))] на 4 сторонах (север, восток, юг, запад) квадрат, который находится внутри. Это означает, что вам никогда не придется искать квадраты по углам. Так как он полностью симметричен c, мы можем взять абсолютные значения входных точек и искать, находится ли точка внутри полуквадратов на положительной стороне координатной плоскости. Это означает, что мы только l oop один раз вместо 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 ((rr / sqrt (2)) * (r / sqrt (2))). Который зацикливается только для половины одного прямоугольника (8-сторонняя симметрия), который находится между внутренним квадратом и внешней границей круга.

...