Skip to content

Instantly share code, notes, and snippets.

@nickponline

nickponline/star_camera.py Secret

Created Jan 23, 2016
Embed
What would you like to do?
# Author: Nicholas Pilkington
# Date: 2016
import cv2
import numpy as np
FOV = 20.0
POV = 1000.0
DEBUG = False
K = 1
class Sky(object):
def __init__(self, x, y, size):
self.x = x
self.y = y
self.size = size
def get_dir(self):
focal = 500.0 / np.tan(10.0 * np.pi / 180.0)
s = [self.x - 499.5, 499.5 - self.y, -focal]
s = np.array(s)
s = s / np.linalg.norm(s)
return s
class Star(object):
def __init__(self, id, ra, dec, mag):
self.id = id
self.ra = ra
self.dec = dec
self.mag = mag
def get_dir(self):
return xyz(self.ra, self.dec)
def angle_between(v1, v2):
v1 = v1 / np.linalg.norm(v1)
v2 = v2 / np.linalg.norm(v2)
dot = np.dot(v1, v2)
dot = min(max(dot, -1.0), 1.0);
return np.arccos(dot)
def xyz(ra, dec):
x = np.cos(ra) * np.cos(dec)
y = np.sin(ra) * np.cos(dec)
z = np.sin(dec)
v = np.array([x, y, z])
v = v / np.linalg.norm(v)
return v
def process_image():
sky = []
# Read image
im = cv2.imread("nightsky.png".format(K), cv2.IMREAD_GRAYSCALE)
ret, im = cv2.threshold(im,127,255,cv2.THRESH_BINARY)
im = 255 - im
params = cv2.SimpleBlobDetector_Params()
# Change thresholds
params.filterByArea = True
params.minArea = 10
# Set up the detector with default parameters.
detector = cv2.SimpleBlobDetector(params)
# Detect blobs.
keypoints = detector.detect(im)
im = 255 - im
for i, kp in enumerate(keypoints):
x, y, size = kp.pt[0], kp.pt[1], kp.size
sk = Sky(x, y, size)
if DEBUG: print x, y, '->', sk.get_dir()
sky.append(sk)
sky = sorted(sky, key=lambda s: -s.size)
for i in xrange(3):
if DEBUG: print i, sky[i].x, sky[i].y, sky[i].size
im_with_keypoints = cv2.drawKeypoints(im, [], np.array([]), (0,255,255), cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
if DEBUG: print len(keypoints)
im_with_keypoints = cv2.resize(im_with_keypoints, (500, 500))
# cv2.imshow("Keypoints", im_with_keypoints)
# cv2.waitKey(0)
return sky
def load_stars():
stars = {}
with open('stars/stars.cat') as fd:
contents = fd.readlines()
for line in contents:
line = line.strip()
line = line.split('\t')
line = map(float, line)
sid, ra, dec, mag = line
stars[int(sid)] = (ra, dec, mag)
return stars
def process_catalogue(threshold=None):
star_dictionary = load_stars()
stars = []
for key, (ra, dec, mag) in star_dictionary.items():
stars.append(Star(key, ra, dec, mag))
stars = sorted(stars, key=lambda star: star.mag)
return stars
def find_transform(indices, stars, blobs):
v1 = stars[indices[0]].get_dir()
v2 = stars[indices[1]].get_dir()
v3 = stars[indices[2]].get_dir()
k1 = blobs[0].get_dir()
k2 = blobs[1].get_dir()
k3 = blobs[2].get_dir()
star_dirs = np.vstack([v1, v2, v3]).T
pic_dirs = np.vstack([k1, k2, k3]).T
inv_star_dirs = np.linalg.inv(star_dirs)
xfm = pic_dirs.dot(inv_star_dirs)
direction =-xfm[2, :]
up = xfm[1, :]
direction = direction / np.linalg.norm(direction)
up = up / np.linalg.norm(up)
return xfm, direction, up
def match(index, indices, stars, sky):
#print 'Stars: ', [stars[k].id for k in indices if k is not None]
if index == len(indices):
xfm, direction, up = find_transform(indices, stars, sky)
det = np.linalg.det(xfm)
if (det < 0.0):
pass
else:
print direction, up
else:
N = len(stars)
for i in xrange(N):
skip = False
if i in indices: continue
for j in xrange(index):
s1 = stars[i]
s2 = stars[indices[j]]
a = angle_between(s1.get_dir(), s2.get_dir())
b = angle_between(sky[j].get_dir(), sky[index].get_dir())
if abs(a - b) > 0.001:
skip = True
if not skip:
indices[index] = i
match(index+1, indices, stars, sky)
indices[index] = None
else:
pass
if __name__ == '__main__':
direction = np.array([1, 0, 0]).T
up = np.array([0, 0, 1]).T
# list of angle, sid, sid smallest star first
stars = process_catalogue(threshold=5)
sky = process_image()
print len(stars), 'stars.'
print len(sky), 'sky objects.'
indices = [None] * 3
match(0, indices, stars, sky)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment