Effectively get lattice points on a squared line
A lattice point is a point with integer coordinates.
The line is the perpendicular divisor between the two grid points A and B. (Each point on this line is equidistant from points A and B.)
How can I efficiently calculate the lattice points on this perpendicular bisector line in the square 0.0 -> N, N?
Here is a square with some examples of points A and B & darr;
Point M is the midpoint between A and B.
My thinking got to me:
Points LA, LB and RA, RB are a square that you can easily calculate to the left and right sides of the AB line.
The LM midpoint between A and LB, and the RM midpoint A and RB are also on the perpendicular bisection line.
So how can you use this information to very quickly compute lattice points on a straight perpendicular line between two points?
this is not homework, simple hobby coding
source to share
I could have thought about this by going matovic in the last draft of the code (which I only took a quick glance at), but anyway ...
Let A = (Ax, Ay), B = (Bx, By), where (Ax, Ay, Bx, By) are integers. Then the line p, perpendicular to the bisector AB, passes through M = (Mx, My) = ((Ax + Bx) / 2, (Ay + By) / 2)
The product of the slopes AB and p is -1, so the slope p is - (Bx - Ax) / (By - Ay)
and, therefore, at the slope, the formula p is (y - My) / (x - Mx) = (Ax - Bx) / (By - Ay)
Reordering,
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
It is clear that for any point of the lattice (x, y) y * (By - Ay) + x * (Bx - Ax) must be an integer. Thus, the line p will pass only through the points of the lattice if (By ^ 2 - Ay ^ 2 + Bx ^ 2 - Ax ^ 2) is even.
Now (By ^ 2 - Ay ^ 2 + Bx ^ 2 - Ax ^ 2) even if and only if (Ax + Bx + Ay + By) is even, which is true if an even number (Ax, Ay, Bx, By ) are odd. In what follows, I will assume that (Ax + Bx + Ay + By) is even.
Let
dx = (Bx - Ax)
dy = (By - Ay)
s = (By ^ 2 - Ay ^ 2 + Bx ^ 2 - Ax ^ 2) / 2
Thus, the equation p is equal to y * dy + x * dx = s
Since y, dy, x, dx, and s are integers, the equation is a linear Diophantine equation, and the standard way to find solutions to such an equation is to use the extended Euclidean algorithm . Our equation will only have solutions if gcd (greatest common divisor) dx and dy divides s. Fortunately, this is true in this case, but I will not provide evidence here.
Let Y, X be the solution y * dy + x * dx = g, where g is 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, therefore Y * dy '+ X * dx' = 1
Dividing the last equation for p by g, we get y * dy '+ x * dx' = s'
And now we can build one solution for it.
(Y * s') * dy '+ (X * s') * dx' = s'
i.e. (X * s ', Y * s') is a lattice point on a line.
We can get all such solutions:
(Y * s' + k * dx ') * dy' + (X * s' - k * dy ') * dx' = s', for all integer k.
To restrict grid solutions 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 will not show solutions to these inequalities right here; see the code below for details.
#! /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 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 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 works correctly. :)
source to share
Okay, I didn't explain my solution exactly, let's start over. Given a doubled resolution grid, the midpoint M will be on the grid. The minimum direction vector of the perpendicular bisector is given by the formula V = (yB - yA, xA - xB) / gcd (yB - yA, xA - xB). Then we look at M and V modulo Z / 2Z x Z / 2Z to see if we can find a point M + iV with even coordinates (for example, on a coarse grid). We can then compute the starting point S = M + jV (j = 0 or 1) on the lattice and get the famous set of points as {S + iV, i is integer}.
[Now running;)] This C ++ code prints S and V as the nearest lattice point to the middle and a vector that can be added or subtracted to get the next / previous lattice point. Then you have to filter the points to get them inside the square (check 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;
}
}
source to share