celutz/panorama/management/commands/test_projections.py

# -*- coding: utf-8 -*-

from __future__ import unicode_literals, division, print_function

from math import log, tan, pi, radians

import matplotlib.pyplot as plt
import numpy as np

from django.core.management.base import BaseCommand, CommandError
from django.conf import settings

from panorama.models import Panorama, Reference

def bearing_diff(b1, b2):
    """In degrees"""
    return (b1 - b2) % 360

def projection_factory(transform):
    def projection(location, p1, p2):
        """For now, simply returns the scaling factors for x and y given two reference points"""
        dx = (p2.x - p1.x) / (radians(bearing_diff(location.bearing(p2.reference_point), location.bearing(p1.reference_point))))
        dy = (p2.y - p1.y) / (transform(location.elevation(p2.reference_point)) - transform(location.elevation(p1.reference_point)))
        return (dx, dy)
    return projection

equirectangular = projection_factory(lambda phi: radians(phi))
cylindrical = projection_factory(lambda phi: tan(radians(phi)))
mercator = projection_factory(lambda phi: log(tan(pi/4 + radians(phi)/2)))
projections = [('equirectangular', equirectangular),
               ('cylindrical', cylindrical),
               ('mercator', mercator)]

class Command(BaseCommand):

    def add_arguments(self, parser):
        parser.add_argument('pano_id', type=int)

    def handle(self, *args, **options):
        p = Panorama.objects.get(pk=options['pano_id'])
        for proj_name, proj in projections:
            self.stdout.write("\n{}".format(proj_name))
            data = {'distx': [], 'disty': [], 'dx': [], 'dy': [], 'ref1': [], 'ref2': [], 'middlex': [], 'ref1.x': []}
            for ref1 in p.panorama_references.order_by('x'):
                for ref2 in p.panorama_references.order_by('x'):
                    if ref2 == ref1:
                        continue
                    if ref2.x >= ref1.x:
                        (dx, dy) = proj(p, ref1, ref2)
                    else:
                        (dx, dy) = proj(p, ref2, ref1)
                    distx = ref2.x - ref1.x
                    disty = ref2.y - ref1.y
                    middlex = (ref2.x + ref1.x) / 2
                    data['ref1'].append(ref1)
                    data['ref2'].append(ref2)
                    data['ref1.x'].append(ref1.x)
                    data['distx'].append(distx)
                    data['disty'].append(disty)
                    data['middlex'].append(middlex)
                    data['dx'].append(dx)
                    data['dy'].append(dy)
                    #self.stdout.write('{} - {}'.format(ref1.reference_point.name, ref2.reference_point.name))
                    #self.stdout.write('dx = {}'.format(dx))
                    #self.stdout.write('dy = {}'.format(dy))
            # Detect outliers
            for axis in ['dx', 'dy']:
                median = np.median(data[axis])
                self.stdout.write("\nMedian for {}: {}".format(axis, median))
                for ref1, ref2, d in zip(data['ref1'], data['ref2'], data[axis]):
                    if ref2.pk < ref1.pk:
                        continue
                    if abs(d - median) / median >= 0.15:
                        self.stdout.write("Outlier: ({:5}, {:5}) - ({:5}, {:5}) → {}={}".format(ref1.x, ref1.y, ref2.x, ref2.y, axis, d))
            mediandx = np.median(data['dx'])
            mediandy = np.median(data['dy'])
            if proj_name == 'equirectangular':
                for xvar in ['distx', 'disty', 'middlex', 'ref1.x']:
                    fig, ax = plt.subplots()
                    ax.scatter(x=data[xvar], y=data['dx'], c=data['ref1.x'], alpha=0.5)
                    ax.hlines([mediandx], 0, 1, transform=ax.get_yaxis_transform(), colors='r')
                    ax.set_title('dx as a function of {}'.format(xvar))
                    plt.show()
            for xvar in ['dx', 'distx', 'disty', 'middlex', 'ref1.x']:
                fig, ax = plt.subplots()
                ax.scatter(x=data[xvar], y=data['dy'], c=data['ref1.x'], alpha=0.5)
                ax.hlines([mediandy], 0, 1, transform=ax.get_yaxis_transform(), colors='r')
                ax.set_title('dy as a function of {}, for {} projection'.format(xvar, proj_name))
                plt.show()