环面上均匀（但不是随机）间隔的 n 点集
Uniformly (but not random) spaced set of n-points on an annulus
美好的一天，
我目前在工作中面临着一个挑战，解决起来非常痛苦。
挑战：
我需要为将被激光切割的穿孔环形板创建孔型的 CAD 绘图。这些孔需要尽可能接近等距，以允许均匀的气流分布。孔的总面积需要等于总环形面积的 40%。根据我的计算，我需要 10,190 个 6.35 毫米的孔来完成这个...
解决方法：
编写一个小程序来计算点集的 XY 坐标，我可以将其导入到我的 CAD 软件中。
问题：
为此，我想使用 Fermat's Spiral approach to calculate the point coordinates. As a starting point, I used info from Marmakoide's Blog。到目前为止，我已经均匀地完成了点分布，并且看起来大致等距分布。我的问题是我需要以某种方式指定 N 点的点集必须落在其中的圆环的内半径和外半径。请注意，我一点也不擅长数学，所以请尽可能保持答案清晰。
这是我的代码：
procedure TForm1.btnCalcClick(Sender: TObject);
var golden_angle, radI, radO, rad, theta, x, y : Double;
i, k, holeQTY, index : integer;
xCoords, yCoords : Array of Double;
begin
holeQTY := StrToInt(edtNumHoles.Text);
golden_angle := Pi * (3 - Sqrt(5));
radI := StrToFloat(edtRadInner.Text);
radO := StrToFloat(edtRadOuter.Text);
SetLength(xCoords, holeQTY);
SetLength(yCoords, holeQTY);
for i := 0 to holeQTY - 1 do
begin
theta := i * golden_angle;
rad := Sqrt(i) / Sqrt(holeQTY);
x := rad * Cos(theta);
y := rad * Sin(theta);
xCoords[i] := x;
yCoords[i] := y;
StatusBar1.Panels[1].Text := IntToStr(i+1) + ' of ' + IntToStr(holeQTY);
StatusBar1.Update;
Application.ProcessMessages;
end;
for k := 0 to holeQTY - 1 do
begin
Chart1.Series[0].AddXY(xCoords[k], yCoords[k], '', clBlack);
end;
end;
提前致谢，
马丁
编辑： 考虑到规范化
在Spreading points页面给出的方法是用单位圈来填充。
您必须添加乘法系数 A
rad := A * Sqrt(i) / Sqrt(N);
索引为 N 的最后一个孔应在 radO 内：
A * Sqrt(N) / Sqrt(N) < radO - HoleRadius
so
A = radO - HoleRadius;
现在开始第一个洞：
A * Sqrt(N - holeQTY + 1) / Sqrt(N) > (radI + HoleRadius)
N - holeQTY + 1 > N * Sqr((radI + HoleRadius) / A)
N * (1 - Sqr((radI + HoleRadius) / A)) > holeQTY - 1
N = Ceil((holeQTY - 1) / (1 - Sqr((radI + HoleRadius) / A))
结果：
您必须使用系数 A 绘制索引 (N - holeQTY + 1)..N 的 holeQTY 个点。示例代码：
var
A, golden_angle, radI, radO, rad, theta, MinGap: Double;
N, x, y, i, holeQTY, holeRadius: Integer;
begin
holeQTY := 120;
radI := 50;
radO := 200;
holeRadius := 4;
MinGap := 3;
golden_angle := Pi * (3 - Sqrt(5));
A := Floor(radO - (holeRadius + MinGap));
N := Ceil((holeQTY - 1) / (1 - Sqr((radI + holeRadius + MinGap) / A)));
Canvas.Ellipse(300 - 200, 300 - 200, 300 + 201, 300 + 201);
Canvas.Ellipse(300 - 50, 300 - 50, 300 + 51, 300 + 51);
for i := (N - holeQTY + 1) to N do begin
theta := i * golden_angle;
rad := A * Sqrt(i) / Sqrt(N);
x := 300 + Round(rad * Cos(theta));
y := 300 + Round(rad * Sin(theta));
Canvas.Ellipse(x - holeRadius, y - holeRadius, x + holeRadius + 1, y + holeRadius + 1);
end;
（MinGap介绍前的图片）
双值代码
holeQTY := 17;
radI := 0.1234;
radO := 0.23456;
holeRadius := 0.00635;
MinGap := 0.0000635;
golden_angle := Pi * (3 - Sqrt(5));
A := radO - (holeRadius + MinGap);
N := Ceil((holeQTY - 1) / (1 - Sqr((radI + holeRadius + MinGap) / A)));
for i := (N - holeQTY + 1) to N do begin
theta := i * golden_angle;
rad := A * Sqrt(i) / Sqrt(N);
x := rad * Cos(theta);
y := rad * Sin(theta);
Memo1.Lines.Add(Format('r:%5.4f x:%f y:%f', [rad, x, y]));
end;
产生输出
r:0.1317 x:0.12 y:0.05
r:0.1397 x:-0.13 y:0.05
r:0.1473 x:0.06 y:-0.13
r:0.1545 x:0.05 y:0.15
r:0.1613 x:-0.14 y:-0.08
r:0.1679 x:0.16 y:-0.04
r:0.1742 x:-0.10 y:0.14
r:0.1804 x:-0.02 y:-0.18
r:0.1863 x:0.14 y:0.12
r:0.1920 x:-0.19 y:0.01
r:0.1976 x:0.14 y:-0.14
r:0.2030 x:-0.01 y:0.20
r:0.2083 x:-0.13 y:-0.16
r:0.2134 x:0.21 y:0.03
r:0.2184 x:-0.18 y:0.12
r:0.2233 x:0.05 y:-0.22
r:0.2281 x:0.11 y:0.20
最终可能的功能与 GUI 分离
type
TPointDouble = record
X, Y: Double;
end;
function CalcRingPoints(const
InnerRadius,
OuterRadius,
HoleRadius,
CoverageRatio //range 0..1
: Double)
: TArray<TPointDouble>;
//doesn't check input validity and possible hole overlaps!
var
ACoeff, golden_angle, rad, theta, MinGap, Area: Double;
N, i, j, holeQTY: Integer;
begin
holeQTY := Round(CoverageRatio * (Sqr(OuterRadius) - Sqr(InnerRadius)) /
Sqr(HoleRadius));
MinGap := 0.1 * HoleRadius;
golden_angle := Pi * (3 - Sqrt(5));
ACoeff := OuterRadius - (HoleRadius + MinGap);
N := Ceil((holeQTY - 1) / (1 - Sqr((InnerRadius + HoleRadius + MinGap) /
ACoeff)));
SetLength(Result, holeQTY);
for i := (N - holeQTY + 1) to N do begin
theta := i * golden_angle;
rad := ACoeff * Sqrt(i / N);
j := i - (N - holeQTY + 1);
Result[j].X := rad * Cos(theta);
Result[j].Y := rad * Sin(theta);
end;
end;
及其用法
var
Pts: TArray<TPointDouble>;
i: Integer;
begin
Pts := CalcRingPoints(1, 2, 0.2, 0.5);
for i := 0 to High(Pts) do
Memo1.Lines.Add(Format('%d r:%5.4f x:%f y:%f',
[i, Hypot(Pts[i].X, Pts[i].Y), Pts[i].X, Pts[i].Y]));
//37 points
让我们看看。
固态盘的公式很明显，对吧？
r[max] = r_coeff * sqrt(max) 这给了我们 r_coeff = r[max]/sqrt[max];
你的问题好像是当我们的平面不是实心圆盘而是大圆盘减去小圆盘时怎么办。
好吧，愚蠢的懒惰方法是更正孔的数量。由于圆盘的平方与半径的平方成正比，并且由于我们希望我们的孔基本上等距分布，所以我们可以希望当我们有足够多的孔时，孔的份额（百分比）与圆盘的平方成正比，对吗？
所以我们有
N = 10 190 = N_outer_disk - N_inner_diskN_outer_disk / N_inner_disk ≈ S_outer / S_inner = (R_outer / R_inner)^2N_inner_disk / N_outer_disk ≈ S_inner / S_outer = (R_inner / R_outer)^2
而且我们必须运行循环I := N_inner_disk + 1 to N_outer_disk计算
r[I] = r_coeff * SqRt(I) angle[I] = I * 137.508 * degrees-to-radians-coefficient
或者也许我们更好地 运行 来自 N_inner_disk - 10% to N_outer_disk + 10% 的循环，并通过 (hole_radius + 10%) + r_hole[I] < R_outer 和 r_hole[I] > R_inner + (hole_radius + 10%) 过滤所有非常大的输出 - 因为我们的计算是近似的，我们会无论如何都要使用过滤，然后我们可以从头开始使用它。
我相信你还需要添加一个变量 - 孔边缘和 inner/outer 圆边缘之间的最小间隙，这样孔就不会 "drlled"边缘部分超出平台。由于缺少上述变量并且为了简单起见，我在上面将间隙设置为孔半径的 10%，这并不完全正确
回到我们的 N-s...
N = 10 190 ≈ N_outer_disk - N_outer_disk*(R_inner / R_outer)^2 = N_outer_disk * (1 - (R_inner / R_outer)^2)
和
N = 10 190 ≈ N_inner_disk*(R_outer / R_inner)^2 - N_inner_disk = N_inner_disk * ((R_outer / R_inner)^2 - 1)
因此对于给定的 N 值，R_outer、R_inner
N_inner_disk ≈ N / ((R_outer / R_inner)^2 - 1)N_outer_disk ≈ N / (1 - (R_inner / R_outer)^2)
现在 运行 从 N_inner+1 到 N_outer 的循环，检查通过上述几何过滤器的孔总数及其累积平方，如果需要，稍微调整系数稍微好一点 N_inner 和 N_outer 直到你接近适合你的分布。
因此 r_coeff = r[max]/sqrt[max] 会变成 r_coeff ≈ R_inner/sqrt[N_inner] ≈ R_outer/sqrt[N_outer]
检查用这两种方法计算的系数是否彼此非常接近将是另一个检查两个 N-s 本身的计算是否内部一致。
美好的一天，
我目前在工作中面临着一个挑战，解决起来非常痛苦。
挑战：
我需要为将被激光切割的穿孔环形板创建孔型的 CAD 绘图。这些孔需要尽可能接近等距，以允许均匀的气流分布。孔的总面积需要等于总环形面积的 40%。根据我的计算，我需要 10,190 个 6.35 毫米的孔来完成这个...
解决方法：
编写一个小程序来计算点集的 XY 坐标，我可以将其导入到我的 CAD 软件中。
问题：
为此，我想使用 Fermat's Spiral approach to calculate the point coordinates. As a starting point, I used info from Marmakoide's Blog。到目前为止，我已经均匀地完成了点分布，并且看起来大致等距分布。我的问题是我需要以某种方式指定 N 点的点集必须落在其中的圆环的内半径和外半径。请注意，我一点也不擅长数学，所以请尽可能保持答案清晰。
这是我的代码：
procedure TForm1.btnCalcClick(Sender: TObject);
var golden_angle, radI, radO, rad, theta, x, y : Double;
i, k, holeQTY, index : integer;
xCoords, yCoords : Array of Double;
begin
holeQTY := StrToInt(edtNumHoles.Text);
golden_angle := Pi * (3 - Sqrt(5));
radI := StrToFloat(edtRadInner.Text);
radO := StrToFloat(edtRadOuter.Text);
SetLength(xCoords, holeQTY);
SetLength(yCoords, holeQTY);
for i := 0 to holeQTY - 1 do
begin
theta := i * golden_angle;
rad := Sqrt(i) / Sqrt(holeQTY);
x := rad * Cos(theta);
y := rad * Sin(theta);
xCoords[i] := x;
yCoords[i] := y;
StatusBar1.Panels[1].Text := IntToStr(i+1) + ' of ' + IntToStr(holeQTY);
StatusBar1.Update;
Application.ProcessMessages;
end;
for k := 0 to holeQTY - 1 do
begin
Chart1.Series[0].AddXY(xCoords[k], yCoords[k], '', clBlack);
end;
end;
提前致谢，
马丁
编辑： 考虑到规范化
在Spreading points页面给出的方法是用单位圈来填充。
您必须添加乘法系数 A
rad := A * Sqrt(i) / Sqrt(N);
索引为 N 的最后一个孔应在 radO 内：
A * Sqrt(N) / Sqrt(N) < radO - HoleRadius
so
A = radO - HoleRadius;
现在开始第一个洞：
A * Sqrt(N - holeQTY + 1) / Sqrt(N) > (radI + HoleRadius)
N - holeQTY + 1 > N * Sqr((radI + HoleRadius) / A)
N * (1 - Sqr((radI + HoleRadius) / A)) > holeQTY - 1
N = Ceil((holeQTY - 1) / (1 - Sqr((radI + HoleRadius) / A))
结果：
您必须使用系数 A 绘制索引 (N - holeQTY + 1)..N 的 holeQTY 个点。示例代码：
var
A, golden_angle, radI, radO, rad, theta, MinGap: Double;
N, x, y, i, holeQTY, holeRadius: Integer;
begin
holeQTY := 120;
radI := 50;
radO := 200;
holeRadius := 4;
MinGap := 3;
golden_angle := Pi * (3 - Sqrt(5));
A := Floor(radO - (holeRadius + MinGap));
N := Ceil((holeQTY - 1) / (1 - Sqr((radI + holeRadius + MinGap) / A)));
Canvas.Ellipse(300 - 200, 300 - 200, 300 + 201, 300 + 201);
Canvas.Ellipse(300 - 50, 300 - 50, 300 + 51, 300 + 51);
for i := (N - holeQTY + 1) to N do begin
theta := i * golden_angle;
rad := A * Sqrt(i) / Sqrt(N);
x := 300 + Round(rad * Cos(theta));
y := 300 + Round(rad * Sin(theta));
Canvas.Ellipse(x - holeRadius, y - holeRadius, x + holeRadius + 1, y + holeRadius + 1);
end;
（MinGap介绍前的图片）
双值代码
holeQTY := 17;
radI := 0.1234;
radO := 0.23456;
holeRadius := 0.00635;
MinGap := 0.0000635;
golden_angle := Pi * (3 - Sqrt(5));
A := radO - (holeRadius + MinGap);
N := Ceil((holeQTY - 1) / (1 - Sqr((radI + holeRadius + MinGap) / A)));
for i := (N - holeQTY + 1) to N do begin
theta := i * golden_angle;
rad := A * Sqrt(i) / Sqrt(N);
x := rad * Cos(theta);
y := rad * Sin(theta);
Memo1.Lines.Add(Format('r:%5.4f x:%f y:%f', [rad, x, y]));
end;
产生输出
r:0.1317 x:0.12 y:0.05
r:0.1397 x:-0.13 y:0.05
r:0.1473 x:0.06 y:-0.13
r:0.1545 x:0.05 y:0.15
r:0.1613 x:-0.14 y:-0.08
r:0.1679 x:0.16 y:-0.04
r:0.1742 x:-0.10 y:0.14
r:0.1804 x:-0.02 y:-0.18
r:0.1863 x:0.14 y:0.12
r:0.1920 x:-0.19 y:0.01
r:0.1976 x:0.14 y:-0.14
r:0.2030 x:-0.01 y:0.20
r:0.2083 x:-0.13 y:-0.16
r:0.2134 x:0.21 y:0.03
r:0.2184 x:-0.18 y:0.12
r:0.2233 x:0.05 y:-0.22
r:0.2281 x:0.11 y:0.20
最终可能的功能与 GUI 分离
type
TPointDouble = record
X, Y: Double;
end;
function CalcRingPoints(const
InnerRadius,
OuterRadius,
HoleRadius,
CoverageRatio //range 0..1
: Double)
: TArray<TPointDouble>;
//doesn't check input validity and possible hole overlaps!
var
ACoeff, golden_angle, rad, theta, MinGap, Area: Double;
N, i, j, holeQTY: Integer;
begin
holeQTY := Round(CoverageRatio * (Sqr(OuterRadius) - Sqr(InnerRadius)) /
Sqr(HoleRadius));
MinGap := 0.1 * HoleRadius;
golden_angle := Pi * (3 - Sqrt(5));
ACoeff := OuterRadius - (HoleRadius + MinGap);
N := Ceil((holeQTY - 1) / (1 - Sqr((InnerRadius + HoleRadius + MinGap) /
ACoeff)));
SetLength(Result, holeQTY);
for i := (N - holeQTY + 1) to N do begin
theta := i * golden_angle;
rad := ACoeff * Sqrt(i / N);
j := i - (N - holeQTY + 1);
Result[j].X := rad * Cos(theta);
Result[j].Y := rad * Sin(theta);
end;
end;
及其用法
var
Pts: TArray<TPointDouble>;
i: Integer;
begin
Pts := CalcRingPoints(1, 2, 0.2, 0.5);
for i := 0 to High(Pts) do
Memo1.Lines.Add(Format('%d r:%5.4f x:%f y:%f',
[i, Hypot(Pts[i].X, Pts[i].Y), Pts[i].X, Pts[i].Y]));
//37 points
让我们看看。
固态盘的公式很明显，对吧？
r[max] = r_coeff * sqrt(max) 这给了我们 r_coeff = r[max]/sqrt[max];
你的问题好像是当我们的平面不是实心圆盘而是大圆盘减去小圆盘时怎么办。
好吧，愚蠢的懒惰方法是更正孔的数量。由于圆盘的平方与半径的平方成正比，并且由于我们希望我们的孔基本上等距分布，所以我们可以希望当我们有足够多的孔时，孔的份额（百分比）与圆盘的平方成正比，对吗？
所以我们有
N = 10 190 = N_outer_disk - N_inner_diskN_outer_disk / N_inner_disk ≈ S_outer / S_inner = (R_outer / R_inner)^2N_inner_disk / N_outer_disk ≈ S_inner / S_outer = (R_inner / R_outer)^2
而且我们必须运行循环I := N_inner_disk + 1 to N_outer_disk计算
r[I] = r_coeff * SqRt(I)angle[I] = I * 137.508 * degrees-to-radians-coefficient
或者也许我们更好地 运行 来自 N_inner_disk - 10% to N_outer_disk + 10% 的循环，并通过 (hole_radius + 10%) + r_hole[I] < R_outer 和 r_hole[I] > R_inner + (hole_radius + 10%) 过滤所有非常大的输出 - 因为我们的计算是近似的，我们会无论如何都要使用过滤，然后我们可以从头开始使用它。
我相信你还需要添加一个变量 - 孔边缘和 inner/outer 圆边缘之间的最小间隙，这样孔就不会 "drlled"边缘部分超出平台。由于缺少上述变量并且为了简单起见，我在上面将间隙设置为孔半径的 10%，这并不完全正确
回到我们的 N-s...
N = 10 190 ≈ N_outer_disk - N_outer_disk*(R_inner / R_outer)^2 = N_outer_disk * (1 - (R_inner / R_outer)^2)
和
N = 10 190 ≈ N_inner_disk*(R_outer / R_inner)^2 - N_inner_disk = N_inner_disk * ((R_outer / R_inner)^2 - 1)
因此对于给定的 N 值，R_outer、R_inner
N_inner_disk ≈ N / ((R_outer / R_inner)^2 - 1)N_outer_disk ≈ N / (1 - (R_inner / R_outer)^2)
现在 运行 从 N_inner+1 到 N_outer 的循环，检查通过上述几何过滤器的孔总数及其累积平方，如果需要，稍微调整系数稍微好一点 N_inner 和 N_outer 直到你接近适合你的分布。
因此 r_coeff = r[max]/sqrt[max] 会变成 r_coeff ≈ R_inner/sqrt[N_inner] ≈ R_outer/sqrt[N_outer]
检查用这两种方法计算的系数是否彼此非常接近将是另一个检查两个 N-s 本身的计算是否内部一致。