Merge branch '30-perform-data-augmentation' into develop

skais/learn/data_augmentation.py

+import random
+from copy import deepcopy
+
+import numpy as np
+from skais.utils.geometry import vectorize_angle, local_ref, transpose_matrix, \
+    quaternion_multiply, angleize_vector, rotation_quaternions, spherical_to_carthesian, carthesian_to_spherical, \
+    rotate_vector_by_quaternion
+
+base_ref = np.array(
+    [
+        [1, 0, 0],
+        [0, 1, 0],
+        [0, 0, 1]
+    ]
+)
+
+
+def shift_positions(aisPositions, axis_lat, axis_long, angle):
+    copy = deepcopy(aisPositions)
+    rotation_quaternion, inverse_rotation_quaternion = rotation_quaternions(axis_lat, axis_long, angle)
+
+    latitudes = np.radians(aisPositions.df['latitude'].to_numpy())
+    longitudes = np.radians(aisPositions.df['longitude'].to_numpy())
+    cogs = np.radians(aisPositions.df['cog'].to_numpy())
+    headings = np.radians(aisPositions.df['heading'].to_numpy())
+
+    for i in range(len(latitudes)):
+        lat = latitudes[i]
+        long = longitudes[i]
+        v_cog = vectorize_angle(cogs[i])
+        v_heading = vectorize_angle(headings[i])
+
+        new_position = rotate_vector_by_quaternion(spherical_to_carthesian(lat, long), rotation_quaternion, inverse_rotation_quaternion)
+        new_lat, new_long = carthesian_to_spherical(new_position)
+
+        origin_ref = local_ref(lat, long)
+        target_ref = local_ref(new_lat, new_long)
+
+        matrix_ob = transpose_matrix(origin_ref, base_ref)
+        matrix_bt = transpose_matrix(base_ref, target_ref)
+
+        # transpose vectors to base referential
+        v_cog = np.matmul(matrix_ob, v_cog)
+        v_heading = np.matmul(matrix_ob, v_heading)
+
+        # rotate vector
+        v_cog = rotate_vector_by_quaternion(v_cog, rotation_quaternion, inverse_rotation_quaternion)
+        v_heading = rotate_vector_by_quaternion(v_heading, rotation_quaternion, inverse_rotation_quaternion)
+
+        # transpose vector back to local referential of the new position
+        v_cog = np.matmul(matrix_bt, v_cog)
+        v_heading = np.matmul(matrix_bt, v_heading)
+
+        latitudes[i] = new_lat
+        longitudes[i] = new_long
+        cogs[i] = angleize_vector(v_cog)
+        headings[i] = angleize_vector(v_heading)
+
+    copy.df['latitude'] = np.degrees(latitudes)
+    copy.df['longitude'] = np.degrees(longitudes)
+    copy.df['cog'] = np.degrees(cogs)
+    copy.df['heading'] = np.degrees(headings)
+
+    return copy

skais/process/data_augmentation/random_axis_rotation_augmentor.py

+import random
+import numpy as np
+from skais.learn.data_augmentation import shift_positions
+
+
+class RandomAxisRotationEngine:
+    def __init__(self, shifts_per_trajectory=1):
+        self.shifts_per_trajectory = shifts_per_trajectory
+
+    def apply(self, trajectories, shuffle=True):
+        new_trajectories = []
+        for trajectory in trajectories:
+            for _ in range(self.shifts_per_trajectory):
+                new_trajectories.append(shift_positions(trajectory, random.uniform(0, 1)*np.pi - np.pi/2,
+                                                        2*(random.uniform(0, 1)*np.pi) - np.pi,
+                                                        2*(random.uniform(0, 1)*np.pi) - np.pi))
+        if not shuffle:
+            return new_trajectories
+        else:
+            return random.shuffle(new_trajectories)

skais/tests/learn/test_data_augmentation.py

+import unittest
+
+from skais.ais.ais_trajectory import AISTrajectory
+from skais.learn.data_augmentation import *
+import pandas as pd
+
+
+class data_augmentation(unittest.TestCase):
+    def test_shift(self):
+        trajectory = AISTrajectory(
+            pd.DataFrame(
+                {
+                    'ts_sec': [i for i in range(10)],
+                    'latitude': [0 for _ in range(5)] + [45 for _ in range(5)],
+                    'longitude': [0 for _ in range(10)],
+                    'cog': [0 for _ in range(10)],
+                    'heading': [90 for _ in range(10)]
+                }
+            )
+        )
+
+        shifted = shift_positions(trajectory, 0, 0, np.pi / 2)
+        self.assertTrue(True)
+
+    def test_shift_2(self):
+        trajectory = AISTrajectory(
+            pd.DataFrame(
+                {
+                    'ts_sec': [i for i in range(10)],
+                    'latitude': [0 for _ in range(10)],
+                    'longitude': [-150, -120, -60, -30, 0, 30, 60, 120, 150, 180],
+                    'cog': [0 for _ in range(10)],
+                    'heading': [90 for _ in range(10)]
+                }
+            )
+        )
+
+        shifted = shift_positions(trajectory, 0, 0, np.pi / 2)
+        self.assertTrue(True)
+
+    def test_shift_3(self):
+        trajectory = AISTrajectory(
+            pd.DataFrame(
+                {
+                    'ts_sec': [0],
+                    'latitude': [-45],
+                    'longitude': [45],
+                    'cog': [45],
+                    'heading': [120]
+                }
+            )
+        )
+
+        shifted = shift_positions(trajectory, np.pi/2, 0, np.pi / 2)
+        self.assertTrue(True)
+

skais/tests/utils/test_geometry.py

 import unittest
 
-from skais.utils.geometry import bresenham
+import numpy as np
+
+from skais.process.geography import PI
+from skais.utils.geometry import bresenham, vectorize_angle, transpose_matrix, vector_rotation, local_ref
 
 
 class TestGeometry(unittest.TestCase):
@@ -18,7 +21,6 @@ class TestGeometry(unittest.TestCase):
 
         self.assertListEqual(result, expected)
 
-
     def test_bresenham_inverted_2(self):
         result = bresenham(16, 4, 3, 9)
         expected = [(3, 9), (4, 9), (5, 8), (6, 8), (7, 7), (8, 7), (9, 7), (10, 6), (11, 6), (12, 6), (13, 5), (14, 5),
@@ -37,6 +39,74 @@ class TestGeometry(unittest.TestCase):
 
         self.assertListEqual(result, expected)
 
+    def test_vectorize_angle_0(self):
+        vectorized_angle = vectorize_angle(0)
+        np.testing.assert_almost_equal(vectorized_angle, np.array([1, 0, 0]))
+
+    def test_vectorize_angle_half_PI(self):
+        vectorized_angle = vectorize_angle(PI / 2)
+        np.testing.assert_almost_equal(vectorized_angle, np.array([0, 1, 0]))
+
+    def test_vectorize_angle_PI(self):
+        vectorized_angle = vectorize_angle(PI)
+        np.testing.assert_almost_equal(vectorized_angle, np.array([-1, 0, 0]))
+
+    def test_transpose_matrix_1(self):
+        origin = np.array(
+            [[1, 0, 0],
+             [0, 1, 0],
+             [0, 0, 1]]
+        )
+        target = np.array(
+            [[1, 0, 0],
+             [0, 1, 0],
+             [0, 0, 1]]
+        )
+        result = transpose_matrix(origin, target)
+        expected = np.array(
+            [[1, 0, 0],
+             [0, 1, 0],
+             [0, 0, 1]]
+        )
+        np.testing.assert_almost_equal(result, expected)
+
+    def test_transpose_matrix_2(self):
+        origin = np.array(
+            [[1, 0, 0],
+             [0, 1, 0],
+             [0, 0, 1]]
+        )
+        target = np.array(
+            [[0, 1, 0],
+             [1, 0, 0],
+             [0, 0, 1]]
+        )
+        result = transpose_matrix(origin, target)
+        expected = np.array(
+            [[0, 1, 0],
+             [1, 0, 0],
+             [0, 0, 1]]
+        )
+        np.testing.assert_almost_equal(result, expected)
+
+    def test_rotate_vector(self):
+        result = vector_rotation(np.array([0, 0, 1]), 0, 0, np.pi / 2)
+
+        np.testing.assert_almost_equal(result, np.array([0, -1, 0]))
+
+    def test_local_ref_0(self):
+        result = local_ref(0, 0)
+
+        np.testing.assert_almost_equal(result, np.array([[0, 0, -1],
+                                                         [0, 1, 0],
+                                                         [1, 0, 0]]))
+
+    def test_local_ref_1(self):
+        result = local_ref(np.pi/2, 0)
+
+        np.testing.assert_almost_equal(result, np.array([[1, 0, 0],
+                                                         [0, 1, 0],
+                                                         [0, 0, 1]]))
 
 if __name__ == '__main__':
     unittest.main()

skais/utils/geography.py

@@ -4,7 +4,7 @@ from numba import jit
 R = 6371000
 
 
-@jit(nopython=True)
+@jit(nopython=True, fastmath=True)
 def great_circle(lat1, lat2, long1, long2):
     x1 = np.radians(lat1)
     y1 = np.radians(long1)

skais/utils/geometry.py

 import numpy as np
 from numba import jit
+import math as m
+
 
 @jit(nopython=True)
 def dist_on_grid(x1, y1, x2, y2):
@@ -65,3 +67,101 @@ def bresenham(x1, y1, x2, y2):
             pixels.append((x, y))
 
     return pixels
+
+@jit(nopython=True, cache=True)
+def vectorize_angle(angle):
+    return np.array([-np.cos(angle), np.sin(angle), 0])
+
+
+@jit(nopython=True, cache=True)
+def angleize_vector(vector):
+    if vector[2] > 0.0001:
+        raise ValueError
+    return np.arctan2(vector[1], -vector[0])
+
+
+@jit(nopython=True, cache=True)
+def local_ref(lat, long):
+    theta = (np.pi / 2) - lat
+    phi = long
+
+    return np.array(
+        [
+            [np.cos(theta) * np.cos(phi), np.cos(theta) * np.sin(phi), -np.sin(theta)],
+            [-np.sin(phi), np.cos(phi), 0],
+            [np.sin(theta) * np.cos(phi), np.sin(theta) * np.sin(phi), np.cos(theta)]
+        ]
+    )
+
+
+# @jit(nopython=True, cache=True)
+def transpose_matrix(origin_ref, target_ref):
+    i1 = origin_ref[0, :]
+    j1 = origin_ref[1, :]
+    k1 = origin_ref[2, :]
+    i2 = target_ref[0, :]
+    j2 = target_ref[1, :]
+    k2 = target_ref[2, :]
+
+    return np.array(
+        [
+            [i1.dot(i2), j1.dot( i2), k1.dot(i2)],
+            [i1.dot(j2), j1.dot(j2), k1.dot(j2)],
+            [i1.dot(k2), j1.dot(k2), k1.dot(k2)]
+        ]
+    )
+
+    # return np.matmul(origin_ref, target_ref)
+
+
+@jit(nopython=True, cache=True)
+def rotation_quaternions(lat, long, angle):
+    theta = (np.pi / 2) - lat
+    phi = long
+    vx = np.cos(phi) * np.sin(theta)
+    vy = np.sin(phi) * np.sin(theta)
+    vz = np.cos(theta)
+
+    return (np.array([np.cos(angle / 2), vx * np.sin(angle / 2), vy * np.sin(angle / 2), vz * np.sin(angle / 2)]),
+            np.array([np.cos(angle / 2), -vx * np.sin(angle / 2), -vy * np.sin(angle / 2), -vz * np.sin(angle / 2)]))
+
+
+@jit(nopython=True, cache=True)
+def quaternion_multiply(quaternion0, quaternion1):
+    w0, x0, y0, z0 = quaternion0
+    w1, x1, y1, z1 = quaternion1
+    return np.array([(w0 * w1) - (x0 * x1) - (y0 * y1) - (z0 * z1),
+                     (w0 * x1) + (x0 * w1) + (y0 * z1) - (z0 * y1),
+                     (w0 * y1) - (x0 * z1) + (y0 * w1) + (z0 * x1),
+                     (w0 * z1) + (x0 * y1) - (y0 * x1) + (z0 * w1)])
+
+
+@jit(nopython=True, cache=True)
+def vector_rotation(vect, lat, long, angle):
+    q, q_conjugate = rotation_quaternions(lat, long, angle)
+    v = np.insert(vect, 0, 0)
+
+    return quaternion_multiply(quaternion_multiply(q, v), q_conjugate)[1:]
+
+
+@jit(nopython=True, cache=True)
+def spherical_to_carthesian(lat, long):
+    theta = (np.pi / 2) - lat
+    phi = long
+    return np.array([np.cos(phi) * np.sin(theta), np.sin(phi) * np.sin(theta), np.cos(theta)])
+
+
+@jit(nopython=True, cache=True)
+def carthesian_to_spherical(vect):
+    x = vect[0]
+    y = vect[1]
+    z = vect[2]
+    XsqPlusYsq = x ** 2 + y ** 2
+    return m.atan2(z, m.sqrt(XsqPlusYsq)), m.atan2(y, x)
+
+
+@jit(nopython=True, cache=True)
+def rotate_vector_by_quaternion(vector, quaternion, inverse_quaternion):
+    vect = np.zeros(4)
+    vect[1:] = vector
+    return quaternion_multiply(quaternion_multiply(quaternion, vect), inverse_quaternion)[1:]