这个算法用来列举圆内整数点,通过推广可以用来列举圆环内的整数点。
不过现在仅是某些特殊情况,没有推广到一般。限制条件:圆心坐标为整数(可以转换为圆心在原点),半径为正整数。
语言为c++,
点的数据结构为:
stuct Vector2
{
int x;
int y;
};
容器为 std::vector<Vector2> points,算法的目标就是把圆内的所有整数点插入points里。
圆是关于坐标轴对称以及 y=x 和 y=-x 对称的图形所以我们仅需考虑第一象限右下部分(圆的八分之一部分)即可。
设半径为r
第一步,加入原点。因为圆心在原点上。
第二步,加入坐标轴上的点。(1, 0)(2,0)...(r, 0)。
第三步,加入y = x, y = -x 上的点 (1, 1)(2, 2)...(⌊r/1.414⌋, ⌊r/1.414⌋) 注意最后一个是向下取整。
第四步,加入其他的点:
最后的这一部分大概是个扇形:
射线 (2, 1)->(3, 1) 和射线 (2, 1)->(3, 2) 与圆相交的区域。
通过观察发现两个事实:
1.这个区域最低点的y坐标为1,最高点的y坐标为 射线 (2, 1)->(3, 2) 与圆相交的点的 y坐标向下取整。
2.因为圆的半径为整数,所以点(r, 0)在圆上。 (r, 1)在圆外。并且有相当多的点(r - 1, j) 是圆内最右边的整数点。
基于这些事实,我们可以把这个区域分为上下两部分:下半部左边界各点坐标为(i + 1, i),右边界格点的x坐标为r-1。上半部为剩下的点。
将x = r - 1带入圆的方程中求出来的y坐标向下取整,即为下半部分的上边界。
也就是说下半部分是个直角梯形:
图中点A为下半部最高点,点B为上半部最高点
上半部分的左边界与下半部分一样,右边界需要求解:将y坐标带入圆方程,求出x坐标并向下取整,即为右边界。
这样知道了左右范围上下范围,就可用循环将各个点加入到数组当中。当然要记得这只是八分之一,还需要将坐标变换:x、y坐标对调,加入负号等。
下面是完整代码:
#pragma once
#include <cassert>
#include <vector>
#include "vector2.h"
class LatticePoint
{
private:
static void pushAxisPoints(int x, std::vector<Vector2<int>> &points)
{
points.push_back({ x, 0 });
points.push_back({ -x, 0 });
points.push_back({ 0, x });
points.push_back({ 0, -x });
}
static void pushDiagonalPoints(int i, std::vector<Vector2<int>> &points)
{
points.push_back({ i, i });
points.push_back({ i, -i });
points.push_back({ -i, i });
points.push_back({ -i, -i });
}
static void pushQuadrantsPoints(int i, int j, std::vector<Vector2<int>> &points)
{
points.push_back({ i, j });
points.push_back({ i, -j });
points.push_back({ -i, j });
points.push_back({ -i, -j });
points.push_back({ j, i });
points.push_back({ j, -i });
points.push_back({ -j, i });
points.push_back({ -j, -i });
}
public:
//get lattice points in circle:(0, 0) radius
static void getFromCircle(int radius, std::vector<Vector2<int>> &points)
{
assert(radius >= 0);
points.push_back({ 0, 0 });
for (int i = 1; i <= radius; i++)
{
pushAxisPoints(i, points);
}
int x1 = floor(radius / sqrt(2)), y1 = x1;
for (int i = 1; i <= x1; i++)
{
pushDiagonalPoints(i, points);
}
if (radius <= 2)return;
x1++;
if (x1 * x1 + y1 * y1 > radius * radius)
{
y1--;
}
int x2 = radius - 1, y2 = floor(sqrt(radius * radius - x2 * x2));
for (int j = 1; j <= y2 - 1; j++)
{
for (int i = j + 1; i <= x2; i++)
{
pushQuadrantsPoints(i, j, points);
}
}
for (int j = y2; j <= y1; j++)
{
int xe = floor(sqrt(radius * radius - j * j));
for (int i = j + 1; i <= xe; i++)
{
pushQuadrantsPoints(i, j, points);
}
}
}
//get lattice points in circle:(0, 0) (smallradius, bigradius]
static void getFromRing(int bigradius, int smallradius, std::vector<Vector2<int>> &points)
{
assert(bigradius == smallradius + 1 && smallradius >= 0);
pushAxisPoints(bigradius, points);
int x1b = floor(bigradius / sqrt(2));
int x1s = floor(smallradius / sqrt(2));
int x1 = x1b, y1 = x1;
if (x1b > x1s)
pushDiagonalPoints(x1, points);
if (bigradius <= 2)return;
x1++;
if (x1 * x1 + y1 * y1 > bigradius * bigradius)
{
y1--;
}
int x2b = bigradius - 1, y2b = floor(sqrt(bigradius * bigradius - x2b * x2b));
int x2s = smallradius - 1, y2s = floor(sqrt(smallradius * smallradius - x2s * x2s));
for (int j = 1; j <= y2s; j++)
{
pushQuadrantsPoints(x2b, j, points);
}
for (int j = y2b > y2s ? y2b : y2b + 1; j <= y1; j++)
{
int xs = floor(sqrt(smallradius * smallradius - j * j));
xs = xs < j ? j : xs;
int xe = floor(sqrt(bigradius * bigradius - j * j));
for (int i = xs + 1; i <= xe; i++)
{
pushQuadrantsPoints(i, j, points);
}
}
}
static int getNumbleofCircle(double radius)
{
assert(radius >= 0);
unsigned int r = (int)radius;
unsigned int sigma = 0;
for (unsigned int i = 1; i <= r; i++)
{
sigma += (int)sqrt(radius * radius - i * i);
}
return 1 + 4 * r + 4 * sigma;
}
static int getNumberofRing(double bigradius, double smallradius)
{
assert(bigradius > smallradius && smallradius >= 0);
return getNumbleofCircle(bigradius) - getNumbleofCircle(smallradius);
}
};
使用方法:
std::vector<Vector2> list;
LatticePoint::getFromCircle(distance, list);
另外还有仅求圆内整数点数量的函数:int LatticPoint::getNumberofCircle(radius);
根据这个方法我还推导出了求圆环内的整数点的方法,函数为 LatticPoint::getFromRing()。不过注意并不完全是通常意义的圆环:这里求出的是到原点的距离为(r, r+1]的所有整数点,注意这里是左开右闭的区间。你可以根据自己的需要来改动。
这里同样有求圆环内整数点数量的函数:int LatticPoint::getNumberofRing(bigradius, smallradius); 这个同样是左开右闭区间内的点,但不一样的地方是大圆半径只需要大于小圆半径就好了。
最后注意,两个求整数点数量的函数的参数可以是小数。而我提出的两个函数,参数需要为整数。