Saltar al contenido

Algoritmo eficiente para obtener puntos en un círculo alrededor de un centro

Nuestro equipo especializado despúes de varios días de investigación y recopilar de información, encontramos la respuesta, nuestro deseo es que te sea útil para tu proyecto.

Solución:

Muy bien, aquí están los puntos de referencia que prometí.

Configuración

Usé el punto de referencia de Google y la tarea consistía en insertar todos los puntos dentro del perímetro del círculo en un std::vector. Hago un punto de referencia para un conjunto de radios y un centro constante:

radii = 10, 20, 50, 100, 200, 500, 1000
center = 100, 500
  • idioma: C ++ 17
  • compilador: msvc 19.24.28316 x64
  • plataforma: windows 10
  • optimización: O2 (optimización completa)
  • subprocesamiento: ejecución de un solo subproceso

Se prueba la exactitud de los resultados de cada algoritmo (en comparación con la salida del algoritmo OP).

Hasta ahora, se comparan los siguientes algoritmos:

  1. Algoritmo de OP enclosing_square.
  2. Mi algoritmo containing_square.
  3. creativecreatorormaybenot algoritmo edge_walking.
  4. Algoritmo de Mandy007 binary_search.

Resultados

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

Código

El código de prueba completo se encuentra a continuación, puede copiarlo, pegarlo y probarlo usted mismo. fill_circle.cpp contiene la implementación de los diferentes algoritmos.

main.cpp

#include 
#include 
#include 

#include 

#include "fill_circle.hpp"

using namespace std::string_literals;

std::unordered_map bench_tests =

    "enclosing_square", enclosing_square,
    "containing_square", containing_square,
    "edge_walking", edge_walking,
    "binary_search", binary_search,
;

std::vector bench_radii = 10, 20, 50, 100, 200, 500, 1000;

void postprocess(std::vector& points)

    std::sort(points.begin(), points.end());
    //points.erase(std::unique(points.begin(), points.end()), points.end());


std::vector prepare(int radius)

    std::vector 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(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 

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& points);

void enclosing_square(point const& center, int radius, std::vector& points);
void containing_square(point const& center, int radius, std::vector& points);
void edge_walking(point const& center, int radius, std::vector& points);
void binary_search(point const& center, int radius, std::vector& 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& 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& 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& 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& 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);
        
    

Esta es una optimización que reduce 1/4 de la dimensión de búsqueda:

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
    
  

o mejor, mejorando el rendimiento la iteración del segundo círculo for evitando el if condicional

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.*/
  

bueno, creo que otra opción es una búsqueda binaria para el límite superior:

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.*/
  

La idea de la búsqueda binaria es encontrar el límite superior de manera óptima, evitando la if condición y cálculos dentro del for ciclo. Para ello, se comprueba cuál es el entero más grande que marca la distancia entre el punto actual y el radio dentro del círculo.

PD: Lo siento mi inglés.

Código

Basado en la idea de @ScottHunter, se me ocurrió el siguiente algoritmo:

#include 

// 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 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;
    

Algoritmo explicado

Este algoritmo realiza una minuto número de cheques. En concreto, solo comprueba en cada fila hasta que se alcanza el primer punto que forma parte del círculo. Además, saltará puntos a la izquierda del punto previamente identificado en la siguiente fila. Además, al usar simetría, solo la mitad de las filas (n/2 + 1/2 ya que comenzamos en 0) están marcadas.

Visualización

Esta es una visualización del algoritmo que creé. El contorno rojo indica el cuadrado que se habría verificado previamente y los píxeles negros indican el círculo real (con el píxel rojo en el medio siendo el centro). El algoritmo verifica los puntos (marcados en azul) y recorre los puntos válidos (marcados en verde).
Como puede ver, el número de píxeles azules al final es un minuto, es decir, solo hay unos pocos puntos en bucle que no forman parte del círculo. Además, observe que solo el primer píxel verde cada vez necesita una verificación, los otros solo se recorren, por lo que aparecen instantáneamente.

Notas

Los ejes podrían invertirse fácilmente, obviamente.

Esto podría optimizarse aprovechando aún más la simetría, es decir, que las filas serán las mismas que las columnas (pasar por todas las filas es lo mismo que pasar por todas las columnas, de izquierda a derecha, de arriba a abajo, viceversa, vise vera) y bajar solo una cuarta parte de las filas desde el centro sería suficiente para determinar exactamente qué puntos van a formar parte del círculo. Sin embargo, siento que el pequeño aumento de rendimiento que esto va a dar no vale la pena el código adicional.
Si alguien quiere codificarlo, proponga una edición de esta respuesta.

Código con comentarios

#include 

// 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 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;
    

Te mostramos las reseñas y valoraciones de los lectores

Si tienes algún recelo y capacidad de prosperar nuestro enunciado puedes escribir una crónica y con gusto lo analizaremos.

¡Haz clic para puntuar esta entrada!
(Votos: 0 Promedio: 0)



Utiliza Nuestro Buscador

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *