如何在矩形的周边均匀分布点
How to evenly distribute points around the perimeter of a rectangle
我正在寻找一种沿着矩形周边的一部分分布点的方法。这些点之间的距离要均匀。
我有一个矩形（通常是正方形）边界，沿该周界有 2 个点（ps 和 pe）标记点的允许范围。这里我用红色标记了允许的范围：
我需要沿着该线段放置 n 个点（通常是 1-3 个）。这些点需要以 d 的距离均匀分布。所以n0..n1和n1..n2等之间的距离应该都是d。边界点也用于分配目的，因此第一个和最后一个点与 ps/pe 之间的距离也应为 d。
这似乎是一项简单的任务，但我很快意识到这种天真的方法在这里行不通。取线段的长度并除以 n+1 不考虑角点。例如：n = 1，点太靠近 pe:
我的数学很生疏（日常工作通常不需要太多），但我尝试了几种不同的方法并且 none 非常有效。我能够使用向量求解 n = 1，方法是找到 ps 和 pe 之间的中点，找到垂直向量，然后将其与线段相交，如下所示。如果 n 是别的东西，或者即使它可以完成，我也不知道如何使这种方法起作用。
最后一点，如果完全均匀分布是不切实际的，那么一个足够好的近似值就可以了。理想情况下，近似值在整个范围内偏离大致相同的量（而不是说，边缘更差）。
毕达哥拉斯来拯救。
考虑 n = 2 有一个角 1 的情况。其余情况可能与更多的边缘情况和可能的配置相似 2.
让我们命名垂直段的长度i和水平段的长度j。
我们要找一个点X，Y，这样距离dist(ps,x) = dist(x,y) = dist(y,pe)
首先，假设角点位于 X 点和 Y 点之间。在这种情况下，我们正在寻找以下方程的解：
x2 = (i - x)2 + (j - x)2 其中 i 和 j 已知。
如果上述方程无解，则意味着角点必须位于ps和X之间或Y和pe之间。我只会介绍 ps 和 X 的情况，因为另一个是对称的：
x2 = i2 + (j-2x)2 是这种情况下的等式。
点和角越多，点和角的可能配置就越多。然而，由于距离应该相等并且长度已知，一系列二次方程很可能就足够了。不过，5+ 分和三个角球的情况有点不确定。
1 n=1 的情况与我将介绍的 n=2 的非对称变体几乎相同。
2这会让它们变得很丑。
任一侧的点根据题目要求等距。只需猜测距离并使用二进制搜索来找到任何精度的解决方案。确定每个矩形边上有多少个点有一个稍微棘手的部分。 （希望只是“轻微”。:)
我打算提出一个算法，但是由于推导在数学上有点乱，我没有足够的时间仔细考虑并仔细检查其正确性。另外，我可能已经包括了一些冗余检查，如果证明某些适当的不等式可能会变得多余，并且可以证明在合理假设下可能始终存在解决方案的存在。我相信这个想法是正确的，但我写这个东西可能会犯一些错误，所以要小心。
因为根据你的评论，在正方形边界线段内只有一个角就足以解决由于对称性导致的其余情况，我将重点关注一个角的情况。
您的一个 90 度角的多边形线段被分成一对垂直的直线段，第一个长度 l1，第二个长度 l2。这两个长度给你。您还想在总长度为 l1 + l2 的多边形线段上添加给定数量的 n 点，以便任意两个连续点之间的欧氏直线距离相同。调用那个未知的距离 d。当你这样做时，你将在 l1 上得到 n1 未知长度的完整片段 d 和 n2 未知长度的完整片段 d l2 这样
n1 + n2 = n
一般来说，您还会在 l1 上得到一段额外的长度 d1 <= d，一端位于 90 度角。类似地，您还将在 l2 上有一段额外的长度 d2 <= d，一端位于 90 度角。因此，d1和d2两条线段共有一端且相互垂直，因此它们形成直角三角形。根据毕达哥拉斯定理，这两段满足方程
d1^2 + d2^2 = d^2
如果我们结合到目前为止的所有方程和信息，我们得到一个方程和限制的系统是：
n1*d + d1 = l1
n2*d + d2 = l2
d1^2 + d2^2 = d^2
n1 + n2 = n
n1 and n2 are non-negative integers
其中变量 d, d1, d2, n1, n2 未知，而 l1, l2, n 已给出。
从前两个等式，您可以表达 d1 和 d2 并代入第三个等式：
d1 = l1 - n1*d
d2 = l2 - n2*d
(l1 - n1*d)^2 + (l2 - n2*d)^2 = d^2
n1 + n2 = n
n1 and n2 are non-negative integers
在特殊情况下，当一个人只想加一个点时，即 n = 1，根据 l1 > l2 或 n1 = n = 1 或 n2 = n = 1 =54=]分别。
说 l1 > l2。然后 n1 = n = 1 和 n2 = 0 所以
d1 = l1 - d
d2 = l2
(l1 - d)^2 + l2^2 = d^2
展开方程，化简求解d：
l1^2 - 2*l1*d + d^2 + l2^2 = d^2
l1^2 + l2^2 - 2*l1*d = 0
d = (l1^2 + l2^2) / (2*l1)
接下来，让我们回到一般情况。你要解决系统
(l1 - n1*d)^2 + (l2 - n2*d)^2 = d^2
n1 + n2 = n
n1 and n2 are non-negative integers
其中变量 d, n1, n2 未知，而 l1, l2, n 已给出。展开第一个等式：
l1^2 - 2 * l1 * n1 * d + n1^2 * d^2 + l2^2 - 2 * l2 * n2 * d + n2^2 * d^2 = d^2
n1 + n2 = n
n1 and n2 are non-negative integers
并将术语组合在一起
(n1^2 + n2^2 - 1) * d^2 - 2 * (l1*n1 + l2*n2) * d + (l1^2 + l2^2) = 0
n1 + n2 = n
n1 and n2 are non-negative integers
第一个方程是d
中的二次方程
(n1^2 + n2^2 - 1) * d^2 - 2 * (l1*n1 + l2*n2) * d + (l1^2 + l2^2) = 0
根据二次公式，您期望有两个解（通常，您选择正数。
如果两者都为正且d < l1和d < l2，你可能有两种解法：
d = ( (l1*n1 + l2*n2) +- sqrt( (l1*n1 + l2*n2)^2 - (l1^2 + l2^2)*(n1^2 + n2^2 - 1) ) ) / (n1^2 + n2^2 - 1)
n1 + n2 = n
n1 and n2 are non-negative integers
现在，如果你能找到合适的n1和n2，你可以使用上面的二次公式计算出必要的d。
对于存在的解决方案，平方根中的表达式必须是非负的，所以你有不等式限制
d = ( (l1*n1 + l2*n2) +- sqrt( (l1*n1 + l2*n2)^2 - (l1^2 + l2^2)*(n1^2 + n2^2 - 1) ) ) / (n1^2 + n2^2 - 1)
(l1*n1 + l2*n2)^2 - (l1^2 + l2^2)*(n1^2 + n2^2 - 1) >= 0
n1 + n2 = n
n1 and n2 are non-negative integers
简化不等式表达：
(l1*n1 + l2*n2)^2 - (l1^2 + l2^2)*(n1^2 + n2^2 - 1) = (l1^2 + l2^2) - (l1*n2 - l2*n1)^2
这将我们带到以下系统
d = ( (l1*n1 + l2*n2) +- sqrt( (l1^2 + l2^2) - (l1*n2 - l2*n1)^2 ) ) / (n1^2 + n2^2 - 1)
(l1^2 + l2^2) - (l1*n2 - l2*n1)^2 >= 0
n1 + n2 = n
n1 and n2 are non-negative integers
分解不平等，
d = ( (l1*n1 + l2*n2) +- sqrt( (l1^2 + l2^2) - (l1*n2 - l2*n1)^2 ) ) / (n1^2 + n2^2 - 1)
(sqrt(l1^2 + l2^2) - l1*n2 + l2*n1) * (sqrt(l1^2 + l2^2) + l1*n2 - l2*n1) >= 0
n1 + n2 = n
n1 and n2 are non-negative integers
所以这些限制有两种情况：
案例 1：
d = ( (l1*n1 + l2*n2) +- sqrt( (l1^2 + l2^2) - (l1*n2 - l2*n1)^2 ) ) / (n1^2 + n2^2 - 1)
sqrt(l1^2 + l2^2) - l1*n2 + l2*n1 >= 0
sqrt(l1^2 + l2^2) + l1*n2 - l2*n1 >= 0
n1 + n2 = n
n1 and n2 are positive integers
或
案例 2：
d = ( (l1*n1 + l2*n2) +- sqrt( (l1^2 + l2^2) - (l1*n2 - l2*n1)^2 ) ) / (n1^2 + n2^2 - 1)
sqrt(l1^2 + l2^2) - l1*n2 + l2*n1 <= 0
sqrt(l1^2 + l2^2) + l1*n2 - l2*n1 <= 0
n1 + n2 = n
n1 and n2 are positive integers
我们关注案例 1，发现案例 2 是不可能的。首先表达 n2 = n - n1，然后将其代入两个不等式中的每一个，并在每个不等式的一侧隔离 n1。此过程产生：
案例 1：
d = ( (l1*n1 + l2*n2) +- sqrt( (l1^2 + l2^2) - (l1*n2 - l2*n1)^2 ) ) / (n1^2 + n2^2 - 1)
( l1*n - sqrt(l1^2 + l2^2) ) / (l1 + l2) <= n1 <= ( l1*n + sqrt(l1^2 + l2^2) ) / (l1 + l2)
n1 + n2 = n
n1 and n2 are positive integers
可以看出情况2颠倒了不等式，这是不可能的，因为左边小于右边
所以算法可能是这样的：
function d = find_d(l1, l2, n)
{
if n = 1 and l1 > l2 {
return d = (l1^2 + l2^2) / (2*l1)
} else if n = 1 and l1 <= l2 {
return d = (l1^2 + l2^2) / (2*l2)
}
for integer n1 >= 0 starting from floor( ( l1*n - sqrt(l1^2 + l2^2) ) / (l1 + l2) ) to floor( ( l1*n + sqrt(l1^2 + l2^2) ) / (l1 + l2) ) + 1
{
n2 = n - n1
D = (l1^2 + l2^2) - (l1*n2 - l2*n1)^2
if D >= 0
{
d1 = ( (l1*n1 + l2*n2) - sqrt( D ) ) / (n1^2 + n2^2 - 1)
d2 = ( (l1*n1 + l2*n2) + sqrt( D ) ) / (n1^2 + n2^2 - 1)
if 0 < d1 < max(l1, l2) {
return d1
} else if 0 < d2 < max(l1, l2) {
return d2
} else {
return "could not find a solution"
}
}
}
}
这是初步版本，建议谨慎使用。我没有足够的时间来检查算法是否可能存在一些接近退化的情况，为此可能必须在某个地方添加一些带有 if 语句的短循环。但是，一般来说，这可能几乎总是有效。我只是发布了一个 python 实现，但是当我找到更多时间并且如果你想让我这样做时，我可以写下这个算法背后的数学原理。这个算法的一些思想可以简化之前的一个角
import numpy as np
import math
import matplotlib.pyplot as plt
def sq_root(x, m, K):
return math.sqrt(x**2 - (K - m*x)**2)
def f(x, n, L):
return sq_root(x, n[0], L[0]) + sq_root(x, n[2], L[2]) + n[1]*x - L[1]
def df(x, n, L):
return ((1-n[0]**2)*x + L[0]*n[0])/sq_root(x, n[0], L[0]) + ((1-n[2]**2)*x + L[2]*n[2])/sq_root(x, n[2], L[2]) + n[1]
#Solving the nonlinear equation for d by using Newton's method:
def solve_f(n, L):
x = sum(L) / (sum(n) + 2)
y = f(x, n, L)
while abs(y) >= 0.0000001:
x = x - y / df(x, n, L)
y = f(x, n, L)
return x - y / df(x, n, L)
def find_n(L, N):
x0 = sum(L) / (N + 1)
# x <= x0
n = np.array([0,0,0])
n[0] = math.floor(L[0]/x0)
n[2] = math.floor(L[2]/x0)
n[1] = N - n[0] - n[2] - 1
return n
def find_d(L, N):
if N==1:
d2 = (L[2]**2 + L[1]**2 - L[0]**2)/(2*L[1])
return math.sqrt(L[0]**2 + d2**2), np.array([0,0,0])
n = find_n(L, N)
return solve_f(n, L), n
def find_the_points(L, N):
d, n = find_d(L, N)
d2 = math.sqrt(d**2 - (L[0]-n[0]*d)**2)
#d3 = math.sqrt(d**2 - (L[2]-n[2]*d)**2)
p = np.zeros((sum(n) + 3, 2))
p[ 0] = np.array([0, L[1]-L[0]])
p[-1] = np.array([L[1], L[1]-L[2]])
e_x = np.array([1,0])
e_y = np.array([0,1])
corner = np.array([0,L[1]])
for i in range(n[0]):
p[i+1] = p[0] + (i+1)*d*e_y
for i in range(n[1]+1):
p[n[0]+i+1] = corner + (d2 + i*d)*e_x
for i in range(n[2]):
p[-(2+i)] = p[-1] + (i+1)*d*e_y
return p, d, n
'''
Test example:
'''
# lengths of the three straight segments along the edges of a square of edge_length L2:
L1 = 5
L2 = 7
L3 = 3
L = np.array([L1, L2, L3])
N = 7
# N = number of points to be added
# If there are two corners then number of segments aligned with edges of square is N - 1
# total number of equidistant segments is N + 1
# n = n[0], n[1], n[2] represents the number of segments aligned with each
# striaght segment from the rectangular polyline along square's boundary
points, d, n = find_the_points(L, N)
print(points)
print(d)
print(n)
plt.figure()
plt.plot(points[:,0], points[:,1])
for j in range(points.shape[0]):
plt.plot(points[j,0], points[j,1], 'ro')
axx = plt.gca()
axx.set_aspect('equal')
plt.show() # if you need...
我正在寻找一种沿着矩形周边的一部分分布点的方法。这些点之间的距离要均匀。
我有一个矩形（通常是正方形）边界，沿该周界有 2 个点（ps 和 pe）标记点的允许范围。这里我用红色标记了允许的范围：
我需要沿着该线段放置 n 个点（通常是 1-3 个）。这些点需要以 d 的距离均匀分布。所以n0..n1和n1..n2等之间的距离应该都是d。边界点也用于分配目的，因此第一个和最后一个点与 ps/pe 之间的距离也应为 d。
这似乎是一项简单的任务，但我很快意识到这种天真的方法在这里行不通。取线段的长度并除以 n+1 不考虑角点。例如：n = 1，点太靠近 pe:
我的数学很生疏（日常工作通常不需要太多），但我尝试了几种不同的方法并且 none 非常有效。我能够使用向量求解 n = 1，方法是找到 ps 和 pe 之间的中点，找到垂直向量，然后将其与线段相交，如下所示。如果 n 是别的东西，或者即使它可以完成，我也不知道如何使这种方法起作用。
最后一点，如果完全均匀分布是不切实际的，那么一个足够好的近似值就可以了。理想情况下，近似值在整个范围内偏离大致相同的量（而不是说，边缘更差）。
毕达哥拉斯来拯救。
考虑 n = 2 有一个角 1 的情况。其余情况可能与更多的边缘情况和可能的配置相似 2.
让我们命名垂直段的长度i和水平段的长度j。
我们要找一个点X，Y，这样距离dist(ps,x) = dist(x,y) = dist(y,pe)
首先，假设角点位于 X 点和 Y 点之间。在这种情况下，我们正在寻找以下方程的解： x2 = (i - x)2 + (j - x)2 其中 i 和 j 已知。
如果上述方程无解，则意味着角点必须位于ps和X之间或Y和pe之间。我只会介绍 ps 和 X 的情况，因为另一个是对称的：
x2 = i2 + (j-2x)2 是这种情况下的等式。
点和角越多，点和角的可能配置就越多。然而，由于距离应该相等并且长度已知，一系列二次方程很可能就足够了。不过，5+ 分和三个角球的情况有点不确定。
1 n=1 的情况与我将介绍的 n=2 的非对称变体几乎相同。
2这会让它们变得很丑。
任一侧的点根据题目要求等距。只需猜测距离并使用二进制搜索来找到任何精度的解决方案。确定每个矩形边上有多少个点有一个稍微棘手的部分。 （希望只是“轻微”。:)
我打算提出一个算法，但是由于推导在数学上有点乱，我没有足够的时间仔细考虑并仔细检查其正确性。另外，我可能已经包括了一些冗余检查，如果证明某些适当的不等式可能会变得多余，并且可以证明在合理假设下可能始终存在解决方案的存在。我相信这个想法是正确的，但我写这个东西可能会犯一些错误，所以要小心。
因为根据你的评论，在正方形边界线段内只有一个角就足以解决由于对称性导致的其余情况，我将重点关注一个角的情况。
您的一个 90 度角的多边形线段被分成一对垂直的直线段，第一个长度 l1，第二个长度 l2。这两个长度给你。您还想在总长度为 l1 + l2 的多边形线段上添加给定数量的 n 点，以便任意两个连续点之间的欧氏直线距离相同。调用那个未知的距离 d。当你这样做时，你将在 l1 上得到 n1 未知长度的完整片段 d 和 n2 未知长度的完整片段 d l2 这样
n1 + n2 = n
一般来说，您还会在 l1 上得到一段额外的长度 d1 <= d，一端位于 90 度角。类似地，您还将在 l2 上有一段额外的长度 d2 <= d，一端位于 90 度角。因此，d1和d2两条线段共有一端且相互垂直，因此它们形成直角三角形。根据毕达哥拉斯定理，这两段满足方程
d1^2 + d2^2 = d^2
如果我们结合到目前为止的所有方程和信息，我们得到一个方程和限制的系统是：
n1*d + d1 = l1
n2*d + d2 = l2
d1^2 + d2^2 = d^2
n1 + n2 = n
n1 and n2 are non-negative integers
其中变量 d, d1, d2, n1, n2 未知，而 l1, l2, n 已给出。
从前两个等式，您可以表达 d1 和 d2 并代入第三个等式：
d1 = l1 - n1*d
d2 = l2 - n2*d
(l1 - n1*d)^2 + (l2 - n2*d)^2 = d^2
n1 + n2 = n
n1 and n2 are non-negative integers
在特殊情况下，当一个人只想加一个点时，即 n = 1，根据 l1 > l2 或 n1 = n = 1 或 n2 = n = 1 =54=]分别。
说 l1 > l2。然后 n1 = n = 1 和 n2 = 0 所以
d1 = l1 - d
d2 = l2
(l1 - d)^2 + l2^2 = d^2
展开方程，化简求解d：
l1^2 - 2*l1*d + d^2 + l2^2 = d^2
l1^2 + l2^2 - 2*l1*d = 0
d = (l1^2 + l2^2) / (2*l1)
接下来，让我们回到一般情况。你要解决系统
(l1 - n1*d)^2 + (l2 - n2*d)^2 = d^2
n1 + n2 = n
n1 and n2 are non-negative integers
其中变量 d, n1, n2 未知，而 l1, l2, n 已给出。展开第一个等式：
l1^2 - 2 * l1 * n1 * d + n1^2 * d^2 + l2^2 - 2 * l2 * n2 * d + n2^2 * d^2 = d^2
n1 + n2 = n
n1 and n2 are non-negative integers
并将术语组合在一起
(n1^2 + n2^2 - 1) * d^2 - 2 * (l1*n1 + l2*n2) * d + (l1^2 + l2^2) = 0
n1 + n2 = n
n1 and n2 are non-negative integers
第一个方程是d
(n1^2 + n2^2 - 1) * d^2 - 2 * (l1*n1 + l2*n2) * d + (l1^2 + l2^2) = 0
根据二次公式，您期望有两个解（通常，您选择正数。
如果两者都为正且d < l1和d < l2，你可能有两种解法：
d = ( (l1*n1 + l2*n2) +- sqrt( (l1*n1 + l2*n2)^2 - (l1^2 + l2^2)*(n1^2 + n2^2 - 1) ) ) / (n1^2 + n2^2 - 1)
n1 + n2 = n
n1 and n2 are non-negative integers
现在，如果你能找到合适的n1和n2，你可以使用上面的二次公式计算出必要的d。
对于存在的解决方案，平方根中的表达式必须是非负的，所以你有不等式限制
d = ( (l1*n1 + l2*n2) +- sqrt( (l1*n1 + l2*n2)^2 - (l1^2 + l2^2)*(n1^2 + n2^2 - 1) ) ) / (n1^2 + n2^2 - 1)
(l1*n1 + l2*n2)^2 - (l1^2 + l2^2)*(n1^2 + n2^2 - 1) >= 0
n1 + n2 = n
n1 and n2 are non-negative integers
简化不等式表达：
(l1*n1 + l2*n2)^2 - (l1^2 + l2^2)*(n1^2 + n2^2 - 1) = (l1^2 + l2^2) - (l1*n2 - l2*n1)^2
这将我们带到以下系统
d = ( (l1*n1 + l2*n2) +- sqrt( (l1^2 + l2^2) - (l1*n2 - l2*n1)^2 ) ) / (n1^2 + n2^2 - 1)
(l1^2 + l2^2) - (l1*n2 - l2*n1)^2 >= 0
n1 + n2 = n
n1 and n2 are non-negative integers
分解不平等，
d = ( (l1*n1 + l2*n2) +- sqrt( (l1^2 + l2^2) - (l1*n2 - l2*n1)^2 ) ) / (n1^2 + n2^2 - 1)
(sqrt(l1^2 + l2^2) - l1*n2 + l2*n1) * (sqrt(l1^2 + l2^2) + l1*n2 - l2*n1) >= 0
n1 + n2 = n
n1 and n2 are non-negative integers
所以这些限制有两种情况：
案例 1：
d = ( (l1*n1 + l2*n2) +- sqrt( (l1^2 + l2^2) - (l1*n2 - l2*n1)^2 ) ) / (n1^2 + n2^2 - 1)
sqrt(l1^2 + l2^2) - l1*n2 + l2*n1 >= 0
sqrt(l1^2 + l2^2) + l1*n2 - l2*n1 >= 0
n1 + n2 = n
n1 and n2 are positive integers
或 案例 2：
d = ( (l1*n1 + l2*n2) +- sqrt( (l1^2 + l2^2) - (l1*n2 - l2*n1)^2 ) ) / (n1^2 + n2^2 - 1)
sqrt(l1^2 + l2^2) - l1*n2 + l2*n1 <= 0
sqrt(l1^2 + l2^2) + l1*n2 - l2*n1 <= 0
n1 + n2 = n
n1 and n2 are positive integers
我们关注案例 1，发现案例 2 是不可能的。首先表达 n2 = n - n1，然后将其代入两个不等式中的每一个，并在每个不等式的一侧隔离 n1。此过程产生：
案例 1：
d = ( (l1*n1 + l2*n2) +- sqrt( (l1^2 + l2^2) - (l1*n2 - l2*n1)^2 ) ) / (n1^2 + n2^2 - 1)
( l1*n - sqrt(l1^2 + l2^2) ) / (l1 + l2) <= n1 <= ( l1*n + sqrt(l1^2 + l2^2) ) / (l1 + l2)
n1 + n2 = n
n1 and n2 are positive integers
可以看出情况2颠倒了不等式，这是不可能的，因为左边小于右边
所以算法可能是这样的：
function d = find_d(l1, l2, n)
{
if n = 1 and l1 > l2 {
return d = (l1^2 + l2^2) / (2*l1)
} else if n = 1 and l1 <= l2 {
return d = (l1^2 + l2^2) / (2*l2)
}
for integer n1 >= 0 starting from floor( ( l1*n - sqrt(l1^2 + l2^2) ) / (l1 + l2) ) to floor( ( l1*n + sqrt(l1^2 + l2^2) ) / (l1 + l2) ) + 1
{
n2 = n - n1
D = (l1^2 + l2^2) - (l1*n2 - l2*n1)^2
if D >= 0
{
d1 = ( (l1*n1 + l2*n2) - sqrt( D ) ) / (n1^2 + n2^2 - 1)
d2 = ( (l1*n1 + l2*n2) + sqrt( D ) ) / (n1^2 + n2^2 - 1)
if 0 < d1 < max(l1, l2) {
return d1
} else if 0 < d2 < max(l1, l2) {
return d2
} else {
return "could not find a solution"
}
}
}
}
这是初步版本，建议谨慎使用。我没有足够的时间来检查算法是否可能存在一些接近退化的情况，为此可能必须在某个地方添加一些带有 if 语句的短循环。但是，一般来说，这可能几乎总是有效。我只是发布了一个 python 实现，但是当我找到更多时间并且如果你想让我这样做时，我可以写下这个算法背后的数学原理。这个算法的一些思想可以简化之前的一个角
import numpy as np
import math
import matplotlib.pyplot as plt
def sq_root(x, m, K):
return math.sqrt(x**2 - (K - m*x)**2)
def f(x, n, L):
return sq_root(x, n[0], L[0]) + sq_root(x, n[2], L[2]) + n[1]*x - L[1]
def df(x, n, L):
return ((1-n[0]**2)*x + L[0]*n[0])/sq_root(x, n[0], L[0]) + ((1-n[2]**2)*x + L[2]*n[2])/sq_root(x, n[2], L[2]) + n[1]
#Solving the nonlinear equation for d by using Newton's method:
def solve_f(n, L):
x = sum(L) / (sum(n) + 2)
y = f(x, n, L)
while abs(y) >= 0.0000001:
x = x - y / df(x, n, L)
y = f(x, n, L)
return x - y / df(x, n, L)
def find_n(L, N):
x0 = sum(L) / (N + 1)
# x <= x0
n = np.array([0,0,0])
n[0] = math.floor(L[0]/x0)
n[2] = math.floor(L[2]/x0)
n[1] = N - n[0] - n[2] - 1
return n
def find_d(L, N):
if N==1:
d2 = (L[2]**2 + L[1]**2 - L[0]**2)/(2*L[1])
return math.sqrt(L[0]**2 + d2**2), np.array([0,0,0])
n = find_n(L, N)
return solve_f(n, L), n
def find_the_points(L, N):
d, n = find_d(L, N)
d2 = math.sqrt(d**2 - (L[0]-n[0]*d)**2)
#d3 = math.sqrt(d**2 - (L[2]-n[2]*d)**2)
p = np.zeros((sum(n) + 3, 2))
p[ 0] = np.array([0, L[1]-L[0]])
p[-1] = np.array([L[1], L[1]-L[2]])
e_x = np.array([1,0])
e_y = np.array([0,1])
corner = np.array([0,L[1]])
for i in range(n[0]):
p[i+1] = p[0] + (i+1)*d*e_y
for i in range(n[1]+1):
p[n[0]+i+1] = corner + (d2 + i*d)*e_x
for i in range(n[2]):
p[-(2+i)] = p[-1] + (i+1)*d*e_y
return p, d, n
'''
Test example:
'''
# lengths of the three straight segments along the edges of a square of edge_length L2:
L1 = 5
L2 = 7
L3 = 3
L = np.array([L1, L2, L3])
N = 7
# N = number of points to be added
# If there are two corners then number of segments aligned with edges of square is N - 1
# total number of equidistant segments is N + 1
# n = n[0], n[1], n[2] represents the number of segments aligned with each
# striaght segment from the rectangular polyline along square's boundary
points, d, n = find_the_points(L, N)
print(points)
print(d)
print(n)
plt.figure()
plt.plot(points[:,0], points[:,1])
for j in range(points.shape[0]):
plt.plot(points[j,0], points[j,1], 'ro')
axx = plt.gca()
axx.set_aspect('equal')
plt.show() # if you need...