nickponline/star_camera.py Secret

## star_camera.py
# 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)
	# 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)