diff --git a/skais/ais/ais_trajectory.py b/skais/ais/ais_trajectory.py index dca3bcc230e78ae872fbcee55bcccc8709a54bdf..68761a7f7a54ea18f874eca3e5e85b39a5622080 100644 --- a/skais/ais/ais_trajectory.py +++ b/skais/ais/ais_trajectory.py @@ -1,5 +1,6 @@ import numbers import random +from copy import copy import pandas as pd import numpy as np @@ -8,7 +9,7 @@ from scipy.interpolate import interp1d from skais.utils.geography import great_circle, position_from_distance, get_coord from skais.ais.ais_points import AISPoints -from skais.utils.geometry import bresenham +from skais.utils.geometry import bresenham, dist_on_grid @jit(nopython=True) @@ -130,22 +131,9 @@ def generate_points(data, positions, height, width, node_size, lower_lat, upper_ data[x, y] = [1] -@jit(nopython=True) -def generate_links(data, positions, height, width, lower_lat, upper_lat, lower_lon, upper_lon): - lon, lat = positions[0, 0], positions[0, 1] - for longitude, latitude in positions[1:]: - x_prv, y_prev = get_coord(lat, lon, height, width, lower_lat, upper_lat, lower_lon, upper_lon) - x_nxt, y_nxt = get_coord(latitude, longitude, height, width, lower_lat, upper_lat, lower_lon, - upper_lon) - - lon, lat = longitude, latitude - for x, y in bresenham(x_prv, y_prev, x_nxt, y_nxt): - data[x, y] = [1] - - @jit(nopython=True) def generate_points_with_features(data, positions, features_vectors, bounds, node_size, height, width, - lower_lat, upper_lat, lower_lon, upper_lon, link): + lower_lat, upper_lat, lower_lon, upper_lon): for pos, f in zip(positions, features_vectors): latitude = pos[1] longitude = pos[0] @@ -162,21 +150,46 @@ def generate_points_with_features(data, positions, features_vectors, bounds, nod for y in range(y_lower_bound, y_upper_bound + 1): for i, v in enumerate(value): data[x, y, i] = v - if link: - lon, lat = positions[0, 0], positions[0, 1] - - value = __get_image_value__(features_vectors[0], bounds) - for pos, f in zip(positions[1:], features_vectors[1:]): - latitude = pos[1] - longitude = pos[0] - x_prv, y_prev = get_coord(lat, lon, height, width, lower_lat, upper_lat, lower_lon, upper_lon) - x_nxt, y_nxt = get_coord(latitude, longitude, height, width, lower_lat, upper_lat, lower_lon, - upper_lon) - lon, lat = longitude, latitude + + +@jit(nopython=True) +def generate_links(data, positions, height, width, lower_lat, upper_lat, lower_lon, upper_lon, values, interpolate=None): + lon, lat = positions[0, 0], positions[0, 1] + + current_value = values[0] + for pos, nxt_value in zip(positions[1:], values[1:]): + latitude = pos[1] + longitude = pos[0] + x_prv, y_prev = get_coord(lat, lon, height, width, lower_lat, upper_lat, lower_lon, upper_lon) + x_nxt, y_nxt = get_coord(latitude, longitude, height, width, lower_lat, upper_lat, lower_lon, + upper_lon) + lon, lat = longitude, latitude + if interpolate is not None and nxt_value != current_value: + dist = dist_on_grid(x_prv, y_prev, x_nxt, y_nxt) for x, y in bresenham(x_prv, y_prev, x_nxt, y_nxt): - for i, v in enumerate(value): - data[x, y, i] = v - value = __get_image_value__(f, bounds) + dist_prev = dist_on_grid(x_prv, y_prev, x, y) + dist_next = dist_on_grid(x, y, x_nxt, y_nxt) + + pixel_color = current_value * (1-dist_prev/dist) + nxt_value * (1-dist_next/dist) + for i in range(len(pixel_color)): + data[x, y, i] = pixel_color[i] + else: + for x, y in bresenham(x_prv, y_prev, x_nxt, y_nxt): + for i in range(len(current_value)): + data[x, y, i] = current_value[i] + current_value = nxt_value + + +@jit(nopython=True) +def generate_values(features_vectors, bounds): + result = np.zeros(features_vectors.shape) + for i in range(len(features_vectors)): + features = features_vectors[i] + value = __get_image_value__(features, bounds) + for j in range(len(value)): + v = value[j] + result[i, j] = v + return result class AISTrajectory(AISPoints): @@ -346,7 +359,7 @@ class AISTrajectory(AISPoints): return result def generate_array_from_positions(self, height=256, width=256, link=True, bounding_box='fit', ref_index=-1, - features=None, node_size=0): + features=None, node_size=0, interpolation=None): positions = self.df[['longitude', 'latitude']].to_numpy() @@ -371,18 +384,23 @@ class AISTrajectory(AISPoints): if features_vectors is not None: nb_channels = len(features_vectors.T) + bounds = np.array(bounds) else: nb_channels = 1 - data = np.zeros((height, width, nb_channels), dtype=np.float) + data = np.zeros((height, width, nb_channels), dtype=float) if features_vectors is None: generate_points(data, positions, height, width, node_size, lower_lat, upper_lat, lower_lon, upper_lon) if link: - generate_links(data, positions, height, width, lower_lat, upper_lat, lower_lon, upper_lon) + generate_links(data, positions, height, width, lower_lat, upper_lat, lower_lon, upper_lon, + np.ones((len(positions), 1))) else: generate_points_with_features(data, positions, features_vectors, np.array(bounds), node_size, height, width, - lower_lat, upper_lat, lower_lon, upper_lon, link) + lower_lat, upper_lat, lower_lon, upper_lon, ) + if link: + generate_links(data, positions, height, width, lower_lat, upper_lat, lower_lon, upper_lon, + generate_values(features_vectors, bounds), interpolation) return data diff --git a/skais/tests/ais/test_ais_trajectory.py b/skais/tests/ais/test_ais_trajectory.py index 51d776ccfe534ef89e1784ad5936dd68b8695572..c491f48c8f0e83a0715b58c6cba655b582677665 100644 --- a/skais/tests/ais/test_ais_trajectory.py +++ b/skais/tests/ais/test_ais_trajectory.py @@ -651,4 +651,32 @@ class TestAISTrajectoryImageGeneration(unittest.TestCase): [[0.25, 1], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0.25, 0.5], [0.5, 0.25]]]) - np.testing.assert_array_equal(result, expected) \ No newline at end of file + np.testing.assert_array_equal(result, expected) + + def test_generate_array_interpolate_links(self): + trajectory = AISTrajectory( + pd.DataFrame( + { + "latitude": [0, 10, 0, 20], + "longitude": [0, 10, 20, 20], + "ts_sec": [i for i in range(4)], + "sog": [10, 10, 20, 40] + } + ) + ) + + result = trajectory.generate_array_from_positions(height=9, width=18, link=True, bounding_box='fit', + features={"sog": (0, 80)}, node_size=0, + interpolation='linear').reshape((9, 18)) + print(result) + expected = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.46875], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.4375], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.40625], + [0, 0, 0, 0, 0, 0, 0, 0.125, 0.125, 0.13526987, 0, 0, 0, 0, 0, 0, 0, 0.375], + [0, 0, 0, 0, 0, 0.125, 0.125, 0, 0, 0, 0.15330406, 0.16458619, 0, 0, 0, 0, 0, 0.34375], + [0, 0, 0, 0.125, 0.125, 0, 0, 0, 0, 0, 0, 0, 0.18154526, 0.19313327, 0, 0, 0, 0.3125], + [0, 0.125, 0.125, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.20959047, 0.22158235, 0, 0.28125], + [0.125, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.23609719, 0.25]]) + + np.testing.assert_array_almost_equal(result, expected) diff --git a/skais/utils/geometry.py b/skais/utils/geometry.py index 7a111d0f32f291f673eff5dca58607c80e8fb446..66f0a29328568a83bd7128674e2f3a9976d29070 100644 --- a/skais/utils/geometry.py +++ b/skais/utils/geometry.py @@ -1,5 +1,12 @@ +import numpy as np from numba import jit +@jit(nopython=True) +def dist_on_grid(x1, y1, x2, y2): + dx = int(x2 - x1) + dy = int(y2 - y1) + return np.sqrt(dx ** 2 + dy ** 2) + @jit(nopython=True) def bresenham(x1, y1, x2, y2): @@ -12,9 +19,8 @@ def bresenham(x1, y1, x2, y2): if dy < 0: sy = -1 - pixels = [] - if abs(dx) > abs(dy): # slope < 1 + if abs(dx) > abs(dy): # slope < 1 if x1 > x2: tmp = x2 x2 = x1 @@ -28,7 +34,7 @@ def bresenham(x1, y1, x2, y2): y = y1 p = (2 * abs(dy)) - abs(dx) pixels.append((x1, y1)) - + for x in range(x1 + 1, x2 + 1): if p < 0: p += 2 * abs(dy) @@ -36,7 +42,7 @@ def bresenham(x1, y1, x2, y2): y += sy p += (2 * abs(dy)) - (2 * abs(dx)) pixels.append((x, y)) - else: # slope >= 1 + else: # slope >= 1 if y1 > y2: tmp = x2 x2 = x1 @@ -59,4 +65,3 @@ def bresenham(x1, y1, x2, y2): pixels.append((x, y)) return pixels -