Lattice point approximation

lattice.py

from math import cos, sin, pi, sqrt, atan
import matplotlib.pyplot as plt
import numpy as np
from mydata import *

# vector addition


def abs2(a):
	return a[0] * a[0] + a[1] * a[1]


def add(a, b):
	return [a[0] + b[0], a[1] + b[1]]


def sub(a, b):
	return [a[0] - b[0], a[1] - b[1]]


def div(a, v):
	return [a[0] / v, a[1] / v]


def mul(a, v):
	return [a[0] * v, a[1] * v]


def rotate(a, theta):
	return [a[0] * cos(theta) - a[1] * sin(theta), a[0] * sin(theta) + a[1] * cos(theta)]


def inside(a):
	return 0 < a[0] < 280 and 0 < a[1] < 280


def flatten(xss):
	return reduce(lambda a, b: a + b, xss)

def gen(a, x, y):
	m = 100
	n = 100
	return filter(inside,
					flatten([[add(a, add(mul(x, i), mul(y, j))) for i in range(-m / 2, m / 2)] for j in
									range(-n / 2, n / 2)])
	)


# Calc lattice from four corner points and number of points on one side.
def make_lattice(a, dist, theta):
	x = rotate([0, dist], theta)
	y = rotate([dist, 0], theta)
	return gen(a, x, y)


def show_points(ps, qs):
	# print(ps)
	ps2 = np.transpose(ps)
	qs2 = np.transpose(qs)
	plt.scatter(ps2[0], ps2[1], color='b', s=5)
	plt.scatter(qs2[0], qs2[1], color='r', s=5)
	plt.axis('equal')
	plt.show(block=True)


def dist2(vs):
	return abs2(sub(vs[0], vs[1]))


def dist(vs):
	return sqrt(dist2(vs))


# ps: theoretical, complete lattice points.
# qs: measured data. May miss some of lattice points.
def calc_diff(ps, qs):
	ds = []
	for q in qs:
		# print ps,q
		r = min([(p, q) for p in ps], key=dist2)
		ds.append(dist2(r))
	return sum(ds)


def angle(a):
	return atan(a[1] / a[0])


def find_initial(qs):
	r = []
	for q in qs:
		qsorted = sorted(qs, key=lambda qq: dist2([qq, q]))[1:5]
		r.append([q, zip(map(lambda a: dist([a, q]), qsorted), qsorted)])
	# print(r)
	ds = reduce(lambda a, b: a + b, map(lambda a: map(lambda b: b[0], a[1]), r))
	ths = reduce(lambda a, b: a + b, map(lambda a: map(lambda b: angle(sub(b[1], a[0])), a[1]), r))
	th = np.mean(filter(lambda a: 0.7 < a < 1.3, ths))
	# plt.hist(ths,np.arange(-pi/2,pi/2,0.02))
	# plt.show()
	# plt.hist(ds,np.arange(10,50,0.25))
	d = np.mean(filter(lambda a: 15 < a < 22, ds))
	# plt.show()
	ps = [qs[len(qs) / 4], qs[len(qs) / 4 * 2], qs[len(qs) / 4 * 3]]
	p = min(ps, key=lambda pp: calc_param_accuracy((pp, d, th)))
	return p, d, th


def calc_param_accuracy(param):
	ps = make_lattice(*param)
	return calc_diff(ps, dat[0])


def perturb(p):
	min_diff = calc_param_accuracy(p)
	x = 0
	y = 0
	for xx in np.arange(-1, 1, 0.2):
		for yy in np.arange(-1, 1, 0.2):
			diff = calc_param_accuracy(([p[0][0] + xx, p[0][1] + yy], p[1], p[2]))
			if diff < min_diff:
				x = xx
				y = yy
				min_diff = diff
				print(diff, x, y)

	return [p[0][0] + x, p[0][1] + y], p[1], p[2]


def main():
	print('start')
	for i in range(1, len(dat)):
		param = find_initial(dat[i])
		print(i, param[0][0], param[0][1], param[1], param[2])

# show_points(make_lattice(*initial_param),dat[i])
#		final = perturb(initial_param)
#		print(initial_param,final)
#		show_points(make_lattice(*final),dat[0])


if __name__ == "__main__":
	main()