2

I've found this answer, which seems to be somewhat related to this question, but I'm wondering if it's possible to generate the coordinates one by one without the additional ~22% (1 - pi / 4) loss of comparing each point to the radius of the circle (by computing the distance between the circle's center and that point).

So far I have the following function in Python. I know by Gauss' circle problem the number of coordinates I will end up with, but I'd like to generate those points one by one as well.

from typing import Iterable
from math import sqrt, floor

def circCoord(sigma: float =1.0, centroid: tuple =(0, 0)) -> Iterable[tuple]:
    r""" Generate all coords within $3\vec{\sigma}$ of the centroid """

    # The number of least iterations is given by Gauss' circle problem:
    # http://mathworld.wolfram.com/GausssCircleProblem.html

    maxiterations = 1 + 4 * floor(3 * sigma) + 4 * sum(\
      floor(sqrt(9 * sigma**2 - i**2)) for i in range(1, floor(3 * sigma) + 1)
    )

    for it in range(maxiterations):
       # `yield` points in image about `centroid` over which we loop

What I'm trying to do is iterate over only those pixels lying within 3 * sigma of a pixel (at centroid in the above function).


I've since written the following example script that demonstrates that the solution below is accurate.

#! /usr/bin/env python3
# -*- coding: utf-8 -*-


import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
import numpy as np
import argparse
from typing import List, Tuple
from math import sqrt


def collect(x: int, y: int, sigma: float =3.0) -> List[Tuple[int, int]]:
    """ create a small collection of points in a neighborhood of some point 
    """
    neighborhood = []

    X = int(sigma)
    for i in range(-X, X + 1):
        Y = int(pow(sigma * sigma - i * i, 1/2))
        for j in range(-Y, Y + 1):
            neighborhood.append((x + i, y + j))

    return neighborhood


def plotter(sigma: float =3.0) -> None:
    """ Plot a binary image """    
    arr = np.zeros([sigma * 2 + 1] * 2)

    points = collect(int(sigma), int(sigma), sigma)

    # flip pixel value if it lies inside (or on) the circle
    for p in points:
        arr[p] = 1

    # plot ellipse on top of boxes to show their centroids lie inside
    circ = Ellipse(\
        xy=(int(sigma), int(sigma)), 
        width=2 * sigma,
        height=2 * sigma,
        angle=0.0
    )

    fig = plt.figure(0)
    ax  = fig.add_subplot(111, aspect='equal')
    ax.add_artist(circ)
    circ.set_clip_box(ax.bbox)
    circ.set_alpha(0.2)
    circ.set_facecolor((1, 1, 1))
    ax.set_xlim(-0.5, 2 * sigma + 0.5)
    ax.set_ylim(-0.5, 2 * sigma + 0.5)

    plt.scatter(*zip(*points), marker='.', color='white')

    # now plot the array that's been created
    plt.imshow(-arr, interpolation='none', cmap='gray')
    #plt.colorbar()

    plt.show()


if __name__ == '__main__':
    parser = argparse.ArgumentParser()

    parser.add_argument('-s', '--sigma', type=int, \
      help='Circle about which to collect points'
    )

    args = parser.parse_args()

    plotter(args.sigma)

And the output for

./circleCheck.py -s 4

is:

enter image description here

Community
  • 1
  • 1
bjd2385
  • 2,013
  • 4
  • 26
  • 47

2 Answers2

5

What about simply something like this (for a circle at origin)?

X = int(R) # R is the radius
for x in range(-X,X+1):
    Y = int((R*R-x*x)**0.5) # bound for y given x
    for y in range(-Y,Y+1):
        yield (x,y)

This can easily be adapted to the general case when the circle is not centred at the origin.

Julien
  • 13,986
  • 5
  • 29
  • 53
0

You might want to consider my algorithm for Gauss circle problem (with some Java source code and an ugly but handy illustration): https://stackoverflow.com/a/42373448/5298879

It's around 3.4x faster than counting points in one of the quarters, plus the center, plus the ones on the axis, that you're doing now, while taking just one more line of code.

You simply imagine an inscribed square and count only one-eighth of what's outside that square inside that circle.

Community
  • 1
  • 1
Zegarek
  • 6,424
  • 1
  • 13
  • 24