I want to get all pixels that are within a circle of a given radius about a given point, where points can only have integer coordinates, i.e. pixels in a canvas.
So I want to obtain all points in the yellow area given (x, y)
and r
.
The most efficient way I can think of is to loop through a square around (x, y)
and check the Euclidean distance for each point:
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.
}
}
}
However, this means that this algorithm will check (r * 2)^2 * (4 - pi) / 4
pixels that are not part of the circle. dx * dx + dy * dy <= r * r
, which seems rather expensive, is called redundantly almost 1 / 4
of the time.
Integrating something like what was proposed here might enhance performance:
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.
}
}
}
However as the author themselves pointed out, this is likely not going to be faster when most of the points are going to be inside of the circle (especially because of abs
), which pi / 4
are in this case.
I was not able to find any resources on this question. I am looking specifically for a solution in C++ and not something in SQL.
Distance Formula for a Point and the Center of a Circle: d=√(x−h)2+(y−k)2 d = ( x − h ) 2 + ( y − k ) 2 , where (x, y) is the point and (h,k) is the center of the circle. This formula is derived from the Pythagorean Theorem.
STEP 1: Take radius as input from the user using std input. STEP 2: Calculate the area of circle using, area = (3.14)*r*r STEP 3: Print the area to the screen using the std output.
For a circle with center (0,0) the equation of a circle is x^2 + y^2 = r^2 (which is the Pythagorean Theorem). So given a center point (h,k) and radius r you can find the points (x,y) that lie on the circumference of the circle.
This lesson will answer this very question. To determine the position of a given point with respect to a circle, all we need to do is to find the distance between the point and the center of the circle, and compare it with the circle's radius. If the distance is greater than the radius, the point lies outside.
Alright here are the benchmarks I promised.
I used google benchmark and the task was to insert all points within the perimiter of the circle into a std::vector<point>
. I benchmark for a set of radii and a constant center:
radii = {10, 20, 50, 100, 200, 500, 1000}
center = {100, 500}
The results of each algorithm is tested for correctness (compared against the output of OPs algorithm).
So far the following algorithms are benchmarked:
enclosing_square
.containing_square
.edge_walking
.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
The complete test code is below, you can copy & paste it and test it yourself. fill_circle.cpp
contains the implementation of the different algorithms.
#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();
}
#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);
#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});
}
}
}
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)
}
}
To avoid repeated generation of pixels at axes, it is worth to start loops from 1 and draw central lines (ix==0 or line==0) in separate loop.
Note that there is also pure integer Bresenham algorithm to generate circumference points.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With