From 61b098c152df19a27159ba66c2ce61dfaf0c44db Mon Sep 17 00:00:00 2001
From: Raphael <raphael.sturgis@gmail.com>
Date: Fri, 5 Aug 2022 14:17:40 +0200
Subject: [PATCH] data augmentation
---
skais/learn/data_augmentation.py | 64 +++++++++++
.../random_axis_rotation_augmentor.py | 20 ++++
skais/tests/learn/__init__.py | 0
skais/tests/learn/test_data_augmentation.py | 56 ++++++++++
skais/tests/utils/test_geometry.py | 74 ++++++++++++-
skais/utils/geography.py | 2 +-
skais/utils/geometry.py | 100 ++++++++++++++++++
7 files changed, 313 insertions(+), 3 deletions(-)
create mode 100644 skais/learn/data_augmentation.py
create mode 100644 skais/process/data_augmentation/random_axis_rotation_augmentor.py
create mode 100644 skais/tests/learn/__init__.py
create mode 100644 skais/tests/learn/test_data_augmentation.py
diff --git a/skais/learn/data_augmentation.py b/skais/learn/data_augmentation.py
new file mode 100644
index 0000000..13da4a4
--- /dev/null
+++ b/skais/learn/data_augmentation.py
@@ -0,0 +1,64 @@
+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
diff --git a/skais/process/data_augmentation/random_axis_rotation_augmentor.py b/skais/process/data_augmentation/random_axis_rotation_augmentor.py
new file mode 100644
index 0000000..f10ebbd
--- /dev/null
+++ b/skais/process/data_augmentation/random_axis_rotation_augmentor.py
@@ -0,0 +1,20 @@
+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)
diff --git a/skais/tests/learn/__init__.py b/skais/tests/learn/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/skais/tests/learn/test_data_augmentation.py b/skais/tests/learn/test_data_augmentation.py
new file mode 100644
index 0000000..3626287
--- /dev/null
+++ b/skais/tests/learn/test_data_augmentation.py
@@ -0,0 +1,56 @@
+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)
+
diff --git a/skais/tests/utils/test_geometry.py b/skais/tests/utils/test_geometry.py
index 7a0b76d..da65b5e 100644
--- a/skais/tests/utils/test_geometry.py
+++ b/skais/tests/utils/test_geometry.py
@@ -1,6 +1,9 @@
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()
diff --git a/skais/utils/geography.py b/skais/utils/geography.py
index ae539e1..99c6e28 100644
--- a/skais/utils/geography.py
+++ b/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)
diff --git a/skais/utils/geometry.py b/skais/utils/geometry.py
index 66f0a29..4178db8 100644
--- a/skais/utils/geometry.py
+++ b/skais/utils/geometry.py
@@ -1,5 +1,7 @@
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:]
--
GitLab