Efficiently get lattice points on line in square

Question

A lattice point is a point with integer coordinates.

The line is the perpendicular bisect between two lattice points A and B. (Every point on that line is equidistant from points A and B.)

How can I efficiently compute the lattice points on that perpendicular bisect line within the square 0,0 → N,N?

Here is a square, with some example points A and B ↓

The point M is the midpoint between A and B.

My thinking has taken me thus far:

The points LA, LB and RA, RB are a square you can easily compute to the left and right sides of the line AB.

The midpoint LM between A and LB, and the midpoint RM A and RB is also on the perpendicular bisect line.

So how can you use this information to very quickly compute the lattice points on the perpendicular bisect line between two points?

_{this isn't homework, its just hobby coding}

Answer 1

I might be over-thinking this, going by matovitch's latest code draft (which I've only had a brief glance at), but anyway...

Let A = (A.x, A.y), B = (B.x, B.y), where (A.x, A.y, B.x, B.y) are integers. Then line p, the perpendicular bisector of AB, passes through
M = (M.x, M.y) = ((A.x + B.x)/2, (A.y + B.y)/2)

The product of the slopes of AB and p is -1, thus the slope of p is
-(B.x - A.x) / (B.y - A.y)
and hence in point-slope form the equation of p is
(y - M.y) / (x - M.x) = (A.x - B.x) / (B.y - A.y)

Rearranging,
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

Clearly, for any lattice point (x, y), y*(B.y - A.y) + x*(B.x - A.x) must be an integer. So the line p will only pass through lattice points if (B.y^2 - A.y^2 + B.x^2 - A.x^2) is even.

Now (B.y^2 - A.y^2 + B.x^2 - A.x^2) is even if & only if (A.x + B.x + A.y + B.y) is even, which is true if an even number of (A.x, A.y, B.x, B.y) are odd. In what follows, I assume that (A.x + B.x + A.y + B.y) is even.

Let
dx = (B.x - A.x)
dy = (B.y - A.y)
s = (B.y^2 - A.y^2 + B.x^2 - A.x^2) / 2
So the equation of p is
y * dy + x * dx = s

Because y, dy, x, dx & s are all integers that equation is a linear Diophantine equation, and the standard way to find the solutions of such an equation is to use the extended Euclidean algorithm. Our equation will only have solutions if the gcd (greatest common divisor) of dx & dy divides s. Fortunately, that's true in this case, but I won't give the proof here.

Let Y, X be a solution of y * dy + x * dx = g, where g is the gcd(dx, dy), i.e.,
Y * dy + X * dx = g
Y * dy/g + X * dx/g = 1

Let dy' = dy/g, dx' = dx/g, s' = s/g, so Y * dy' + X * dx' = 1

Dividing the last equation for p through by g, we get y * dy' + x * dx' = s'

And we can now construct one solution for it.
(Y * s') * dy' + (X * s') * dx' = s'
i.e., (X * s', Y * s') is a lattice point on the line.

We can get all solutions like this:
(Y * s' + k * dx') * dy' + (X * s' - k * dy') * dx' = s', for all integers k.

To restrict the solutions to the grid from (0, 0) to (W, H), we need to solve these inequalities for k:
0 <= X * s' - k * dy' <= W and 0 <= Y * s' + k * dx' <= H

I won't show the solutions of those inequalities right here; for the details see the code below.

#! /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()

I haven't tested this code extensively, but it appears to work correctly. :)

Answer 2

Ok I sure did not explained my solution clearly, let's start again. Given a grid with twice the resolution, the middle point M will be on the grid. The minimal direction vector of the perpendicular bissector is given by V = (yB - yA, xA - xB) / gcd(yB - yA, xA - xB). Then we look at M and V modulo the lattice Z/2Z x Z/2Z to check if one can find a point M + iV with even coordinates (aka on the coarse grid). We can then compute a starting point S = M + jV (j = 0 or 1 in fact) on the lattice and get the famous set of points as {S + iV, i integer}.

[Now running ;)] This C++ code print S and V, aka the nearest lattice point to the middle and the vector one can add or subtract to get the next/previous lattice point. You then have to filter the points to get those inside the square (test it here : 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;
    }
}