gpt4 book ai didi

python - 高效获取正方形线上的格点

转载 作者:太空狗 更新时间:2023-10-30 00:55:06 25 4
gpt4 key购买 nike

格点是具有整数坐标的点。

该线是两个格点 A 和 B 之间的垂直平分线。(该线上的每个点与点 A 和 B 的距离相等。)

如何有效地计算正方形 0,0 → N,N 内垂直平分线上的格点?

这是一个正方形,有一些示例点 A 和 B ↓

enter image description here

点 M 是 A 和 B 之间的中点。

到目前为止,我的思路是:

点 LA、LB 和 RA、RB 是一个正方形,您可以轻松地计算到线 AB 的左侧和右侧。

A和LB的中点LM,A和RB的中点RM也在垂直平分线上。

那么如何使用此信息快速计算两点之间垂直平分线上的格点?

这不是家庭作业,只是爱好编码

最佳答案

根据 matovitch 的最新代码草稿(我只粗略看了一眼),我可能想多了,但无论如何...

令 A = (A.x, A.y), B = (B.x, B.y),其中 (A.x, A.y, B.x, B.y) 是整数。则AB的垂直平分线p穿过
M = (M.x, M.y) = ((A.x + B.x)/2, (A.y + B.y)/2)

AB和p的斜率的乘积是-1,因此p的斜率是
-(B.x - A.x)/(B.y - A.y)
因此,在点斜率形式下,p 的方程是
(y - M.y)/(x - M.x) = (A.x - B.x)/(B.y - A.y)

重新排列,
y*(B.y - A.y) + x*(B.x - A.x) = M.y * (B.y - A.y) + M.x * (B.x - A.x)
= ((B.y + A.y) * (B.y - A.y) + (B.x + A.x) * (B.x - A.x))/2
= (B.y^2 - A.y^2 + B.x^2 - A.x^2)/2

显然,对于任何格点(x, y),y*(B.y - A.y) + x*(B.x - A.x) 必须是一个整数。所以只有当 (B.y^2 - A.y^2 + B.x^2 - A.x^2) 是偶数时,线 p 才会通过格点。

现在 (B.y^2 - A.y^2 + B.x^2 - A.x^2) 是偶数当且仅当 (A.x + B.x + A.y + B.y) 是偶数,如果 (A.x, A.y, B.x, B.y) 是奇数。在下文中,我假设 (A.x + B.x + A.y + B.y) 是偶数。


dx = (B.x - A.x)
dy = (B.y - A.y)
s = (B.y^2 - A.y^2 + B.x^2 - A.x^2)/2
所以p的方程是
y * dy + x * dx = s

因为 y、dy、x、dx 和 s 都是整数,所以方程是线性丢番图方程,求此类方程解的标准方法是使用 extended Euclidean algorithm .只有当 dx 和 dy 的 gcd(最大公约数)除以 s 时,我们的方程才会有解。幸运的是,在这种情况下确实如此,但我不会在这里给出证明。

令 Y, X 为 y * dy + x * dx = g 的解,其中 g 为 gcd(dx, dy),即
Y * dy + X * dx = g
Y * dy/g + X * dx/g = 1

令 dy' = dy/g, dx' = dx/g, s' = s/g, 所以Y * dy' + X * dx' = 1

p 的最后一个方程除以 g,我们得到y * dy' + x * dx' = s'

现在我们可以为它构建一个解决方案。
(Y * s') * dy' + (X * s') * dx' = s'
即,(X * s', Y * s') 是直线上的格点。

我们可以得到这样的所有解决方案:
(Y * s' + k * dx') * dy' + (X * s' - k * dy') * dx' = s',对于所有整数 k。

为了限制从 (0, 0) 到 (W, H) 的网格的解,我们需要解决 k 的这些不等式:
0 <= X * s' - k * dy' <= W 和 0 <= Y * s' + k * dx' <= H

我不会在这里展示这些不等式的解决方案;有关详细信息,请参阅下面的代码。

#! /usr/bin/env python

''' Lattice Line

Find lattice points, i.e, points with integer co-ordinates,
on the line that is the perpendicular bisector of the line segment AB,
where A & B are lattice points.

See http://stackoverflow.com/q/31265139/4014959

Written by PM 2Ring 2015.07.08
Code for Euclid's algorithm & the Diophantine solver written 2010.11.27
'''

from __future__ import division
import sys
from math import floor, ceil

class Point(object):
''' A simple 2D point '''
def __init__(self, x, y):
self.x, self.y = x, y

def __repr__(self):
return "Point(%s, %s)" % (self.x, self.y)

def __str__(self):
return "(%s, %s)" % (self.x, self.y)


def euclid(a, b):
''' Euclid's extended algorithm for the GCD.
Returns a list of tuples of (dividend, quotient, divisor, remainder)
'''
if a < b:
a, b = b, a

k = []
while True:
q, r = a // b, a % b
k.append((a, q, b, r))
if r == 0:
break
a, b = b, r
return k


def dio(aa, bb):
''' Linear Diophantine solver
Returns [x, aa, y, bb, d]: x*aa + y*bb = d
'''
a, b = abs(aa), abs(bb)
swap = a < b
if swap:
a, b = b, a

#Handle trivial cases
if a == b:
eqn = [2, a, -1, a]
elif a % b == 0:
q = a // b
eqn = [1, a, 1-q, b]
else:
#Generate quotients & remainders list
z = euclid(a, b)[::-1]

#Build equation from quotients & remainders
eqn = [0, 0, 1, 0]
for v in z[1:]:
eqn = [eqn[2], v[0], eqn[0] - eqn[2]*v[1], v[2]]

#Rearrange & fix signs, if required
if swap:
eqn = eqn[2:] + eqn[:2]

if aa < 0:
eqn[:2] = [-eqn[0], -eqn[1]]
if bb < 0:
eqn[2:] = [-eqn[2], -eqn[3]]

d = eqn[0]*eqn[1] + eqn[2]*eqn[3]
if d < 0:
eqn[0], eqn[2], d = -eqn[0], -eqn[2], -d

return eqn + [d]


def lattice_line(pA, pB, pC):
''' Find lattice points, i.e, points with integer co-ordinates, on
the line that is the perpendicular bisector of the line segment AB,
Only look for points in the rectangle from (0,0) to C

Let M be the midpoint of AB. Then M = ((A.x + B.x)/2, (A.y + B.y)/2),
and the equation of the perpendicular bisector of AB is
(y - M.y) / (x - M.x) = (A.x - B.x) / (B.y - A.y)
'''

nosolutions = 'No solutions found'

dx = pB.x - pA.x
dy = pB.y - pA.y

#Test parity of co-ords to see if there are solutions
if (dx + dy) % 2 == 1:
print nosolutions
return

#Handle horizontal & vertical lines
if dx == 0:
#AB is vertical, so bisector is horizontal
y = pB.y + pA.y
if dy == 0 or y % 2 == 1:
print nosolutions
return
y //= 2
for x in xrange(pC.x + 1):
print Point(x, y)
return

if dy == 0:
#AB is horizontal, so bisector is vertical
x = pB.x + pA.x
if x % 2 == 1:
print nosolutions
return
x //= 2
for y in xrange(pC.y + 1):
print Point(x, y)
return

#Compute s = ((pB.x + pA.x)*dx + (pB.y + pA.y)*dy) / 2
#s will always be an integer since (dx + dy) is even
#The desired line is y*dy + x*dx = s
s = (pB.x**2 - pA.x**2 + pB.y**2 - pA.y**2) // 2

#Find ex, ey, g: ex * dx + ey * dy = g, where g is the gcd of (dx, dy)
#Note that g also divides s
eqn = dio(dx, dy)
ex, ey, g = eqn[::2]

#Divide the parameters of the equation by the gcd
dx //= g
dy //= g
s //= g

#Find lattice limits
xlo = (ex * s - pC.x) / dy
xhi = ex * s / dy
if dy < 0:
xlo, xhi = xhi, xlo

ylo = -ey * s / dx
yhi = (pC.y - ey * s) / dx
if dx < 0:
ylo, yhi = yhi, ylo

klo = int(ceil(max(xlo, ylo)))
khi = int(floor(min(xhi, yhi)))
print 'Points'
for k in xrange(klo, khi + 1):
x = ex * s - dy * k
y = ey * s + dx * k
assert x*dx + y*dy == s
print Point(x, y)


def main():
if len(sys.argv) != 7:
print ''' Find lattice points, i.e, points with integer co-ordinates,
on the line that is the perpendicular bisector of the line segment AB,
where A & B are lattice points with co-ords (xA, yA) & (xB, yB).
Only print lattice points in the rectangle from (0, 0) to (W, H)

Usage:
%s xA yA xB yB W H''' % sys.argv[0]
exit(1)

coords = [int(s) for s in sys.argv[1:]]
pA = Point(*coords[0:2])
pB = Point(*coords[2:4])
pC = Point(*coords[4:6])
lattice_line(pA, pB, pC)


if __name__ == '__main__':
main()

我还没有广泛测试这段代码,但它似乎可以正常工作。 :)

关于python - 高效获取正方形线上的格点,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/31265139/

25 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com