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)