有效地在正方形上获得格点

Question

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

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

如何有效地计算方形0,0→N,N中垂直平分线上的晶格点？

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

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

到目前为止,我的想法一直带着我:

点LA,LB和RA,RB是可以容易地计算到线AB的左侧和右侧的正方形.

A和LB之间的中点LM,以及中点RM A和RB也在垂直平分线上.

那么如何使用这些信息来快速计算两点之间垂直平分线上的晶格点？

_{这不是功课,它只是爱好编码}

Answer 1

我可能会过度思考这个问题,参加matovitch的最新代码草案(我只是简单地看一眼),但无论如何......

设A =(Ax,Ay),B =(Bx,By),其中(Ax,Ay,Bx,By)是整数.然后线p,AB的垂直平分线,通过
M =(MX,MY)=((AX + Bx的)/ 2,(AY +通过)/ 2)

AB和p的斜率的乘积是-1,因此p的斜率是
- (Bx-Ax)/(By-Ay)
,因此在点斜率形式中p的方程是
(y-My)/( x - Mx)=(Ax - Bx)/(By - Ay)

重新排列,
y*(By-Ay)+ x*(Bx-Ax)= My*(By-Ay)+ Mx*(Bx-Ax)
=((By + Ay)*(By-Ay)+(Bx + Ax)*(Bx-Ax))/ 2
=(By ^ 2-Ay ^ 2 + Bx ^ 2-Ax ^ 2)/ 2

显然,对于任何格点(x,y),y*(By-Ay)+ x*(Bx-Ax)必须是整数.因此,如果(By ^ 2-Ay ^ 2 + Bx ^ 2-Ax ^ 2)是偶数,则线p将仅通过格点.

现在(按^ 2 - Ay ^ 2 + Bx ^ 2 - Ax ^ 2)即使(仅当Ax + Bx + Ay + By)是偶数,如果是偶数(Ax,Ay,Bx)也是如此,By)很奇怪.在下文中,我假设(Ax + Bx + Ay + By)是偶数.

令
dx =(Bx-Ax)
dy =(By-Ay)
s =(By ^ 2-Ay ^ 2 + Bx ^ 2-Ax ^ 2)/ 2
因此p的方程是
y*dy + x*dx =小号

因为y,dy,x,dx和s都是整数,所以方程是线性丢番图方程,找到这种方程的解的标准方法是使用扩展的欧几里德算法.如果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')是该行上的格点.

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

为了将网格的解决方案从(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()

我没有广泛测试此代码,但它似乎正常工作.:)

Answer 2

好吧，我确实没有清楚地解释我的解决方案，让我们重新开始。给定一个两倍分辨率的网格，中点 M 将位于网格上。垂直平分线的最小方向向量由 V = (yB - yA, xA - xB) / gcd(yB - yA, xA - xB) 给出。然后我们查看 M 和 V 对晶格 Z/2Z x Z/2Z 进行模，检查是否可以找到坐标为偶数的点 M + iV（也称为粗网格上）。然后，我们可以计算晶格上的起始点 S = M + jV（实际上 j = 0 或 1），并得到著名的点集 {S + iV, i 整数}。

[正在运行;)] 此 C++ 代码打印 S 和 V，即距离中间最近的格点，以及可以添加或减去以获得下一个/上一个格点的向量。然后，您必须过滤这些点以获取正方形内的点（在这里测试： http: //coliru.stacked-crooked.com/a/ba9f8aec45e1c2ea）：

int gcd(int n1, int n2)
{
    n1 = (n1 > 0) ? n1 : -n1;
    n2 = (n2 > 0) ? n2 : -n2;

    if (n1 > n2) 
    {
        int t = n1;
        n1 = n2;
        n2 = t; 
    } 

    while (n2 % n1 != 0)
    {
        int tmp = n2 % n1;
        n2 = n1;
        n1 = tmp; 
    }

    return n1;
}

struct Point
{

    const Point& operator=(const Point& rhs)
    {
        x = rhs.x;
        y = rhs.y;

        return *this;
    }

    const Point& operator+=(const Point& rhs)
    {
        x += rhs.x;
        y += rhs.y;

        return *this;
    }

    const Point& operator-=(const Point& rhs)
    {
        x += rhs.x;
        y += rhs.y;

        return *this;
    }

    const Point& operator/=(int rhs)
    {
        x /= rhs;
        y /= rhs;

        return *this;
    }

    const Point& reduce()
    {
        return *this /= gcd(x, y);
    }

    int x;
    int y;
};

const Point operator+(Point lhs, const Point& rhs)
{
    lhs += rhs;
    return lhs;
}

const Point operator-(Point lhs, const Point& rhs)
{
    lhs -= rhs;
    return lhs;
}

const Point operator/(Point lhs, int rhs)
{
    lhs /= rhs;
    return lhs;
}

bool findBase(Point& p1, Point& p2, Point& oBase, Point& oDir)
{
    Point m = p1 + p2;
    Point v = {p2.y - p1.y, p1.x - p2.x};

    oDir = v.reduce();

    if (m.x % 2 == 0 && m.y % 2 == 0)
    {
        oBase = m / 2;
        return true;
    } 
    else if (((m.x % 2 == 0 && v.x % 2 == 0) &&
              (m.y % 2 == 1 && v.y % 2 == 1)) ||
             ((m.x % 2 == 1 && v.x % 2 == 1) &&
              (m.y % 2 == 0 && v.y % 2 == 0)) ||
             ((m.x % 2 == 1 && v.x % 2 == 1) &&
              (m.y % 2 == 1 && v.y % 2 == 1)))
    {
        oBase = (m + v) / 2;
        return true;
    }
    else 
    {
        return false;
    }
}