diff --git a/skais/ais/ais_trajectory.py b/skais/ais/ais_trajectory.py
index 9ce6082197096898c856ac87f27287fbb723ca7e..e123aff7898a03ad1eff6d03c94f70b3734d1e71 100644
--- a/skais/ais/ais_trajectory.py
+++ b/skais/ais/ais_trajectory.py
@@ -7,6 +7,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
 
 
 @jit(nopython=True)
@@ -224,7 +225,8 @@ class AISTrajectory(AISPoints):
                 result.append((row['ts_sec'], current_label))
         return result
 
-    def generate_array_from_positions(self, height=256, width=256, link=True, bounding_box='fit', features=None, node_size=0):
+    def generate_array_from_positions(self, height=256, width=256, link=True, bounding_box='fit', features=None,
+                                      node_size=0):
         nb_channels = 1
         if features is not None:
             nb_channels = len(features)
@@ -248,5 +250,12 @@ class AISTrajectory(AISPoints):
                     data[x, y] = [1]
 
         if link:
-            raise ValueError("feature not implemented")
+            lon, lat = positions[0, 0], positions[0, 1]
+            for longitude, latitude in positions[1:]:
+                x_prv, y_prev = get_coord(lat, lon, height, width, min_lat, max_lat, min_lon, max_lon)
+                x_nxt, y_nxt = get_coord(latitude, longitude, height, width, min_lat, max_lat, min_lon, max_lon)
+
+                lon, lat = longitude, latitude
+                for x, y in bresenham(x_prv, y_prev, x_nxt, y_nxt):
+                    data[x, y] = [1]
         return data
diff --git a/skais/tests/ais/test_ais_trajectory.py b/skais/tests/ais/test_ais_trajectory.py
index e4b6c6909161ff8852df6584ca0ed4a2455eb3b8..fe8efc4c087776701515b8036ae9ecf4939e3b27 100644
--- a/skais/tests/ais/test_ais_trajectory.py
+++ b/skais/tests/ais/test_ais_trajectory.py
@@ -428,4 +428,29 @@ class TestAISTrajectory(unittest.TestCase):
                              [1, 1, 0, 0, 0, 0, 0, 0, 0],
                              [1, 1, 0, 0, 0, 0, 0, 0, 0]])
 
+        np.testing.assert_array_equal(result, expected)
+
+    def test_generate_array_from_positions_with_line(self):
+        trajectory = AISTrajectory(
+            pd.DataFrame(
+                {
+                    "latitude": [0, 10, 0, 20],
+                    "longitude": [0, 10, 20, 20],
+                    "ts_sec": [i for i in range(4)]
+                }
+            )
+        )
+
+        result = trajectory.generate_array_from_positions(height=9, width=18, link=True, bounding_box='fit',
+                                                          features=None, node_size=0).reshape((9, 18))
+        expected = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
+                             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
+                             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
+                             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
+                             [0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1],
+                             [0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1],
+                             [0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1],
+                             [0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1],
+                             [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1]])
+
         np.testing.assert_array_equal(result, expected)
\ No newline at end of file
diff --git a/skais/utils/geometry.py b/skais/utils/geometry.py
index 885aa429e8d6b5a5228888da3449d4a8b2ce9f0b..adbd8a58b0d66373ebb2eecd6ca8f3bc8aaecefe 100644
--- a/skais/utils/geometry.py
+++ b/skais/utils/geometry.py
@@ -1,14 +1,4 @@
 def bresenham(x1, y1, x2, y2):
-    if x1 > x2:
-        tmp = x2
-        x2 = x1
-        x1 = tmp
-        
-
-        tmp = y2
-        y2 = y1
-        y1 = tmp
-    
     dx = int(x2 - x1)
     dy = int(y2 - y1)
 
@@ -23,6 +13,14 @@ def bresenham(x1, y1, x2, y2):
 
     pixels = [(x1, y1)]
     if abs(dx) > abs(dy): # slope < 1
+        if x1 > x2:
+            tmp = x2
+            x2 = x1
+            x1 = tmp
+
+            tmp = y2
+            y2 = y1
+            y1 = tmp
         p = (2 * abs(dy)) - abs(dx)
         
         for x in range(x1 + 1, x2 + 1):
@@ -33,6 +31,14 @@ def bresenham(x1, y1, x2, y2):
                 p += (2 * abs(dy)) - (2 * abs(dx))
             pixels.append((x, y))
     else: # slope >= 1
+        if y1 > y2:
+            tmp = x2
+            x2 = x1
+            x1 = tmp
+
+            tmp = y2
+            y2 = y1
+            y1 = tmp
         p = (2 * abs(dx)) - abs(dy)
         for y in range(y1 + 1, y2 + 1):
             if p < 0: