Merge branch 'road'

Enums.py

@@ -0,0 +1,31 @@
+from enum import Enum
+
+
+class DIRECTION(Enum):
+    WEST = 0
+    EAST = 1
+    NORTH = 2
+    SOUTH = 3
+
+
+class COLLUMN_STYLE(Enum):
+    INNER = 1
+    OUTER = 2
+    BOTH = 3
+
+
+class LINE_OVERLAP(Enum):
+    NONE = 0
+    MAJOR = 1
+    MINOR = 2
+
+
+class LINE_THICKNESS_MODE(Enum):
+    MIDDLE = 0
+    DRAW_COUNTERCLOCKWISE = 1
+    DRAW_CLOCKWISE = 2
+
+
+class ROTATION(Enum):
+    CLOCKWISE = 0
+    COUNTERCLOCKWISE = 1

main.py

@@ -5,11 +5,11 @@ import gdpc.exceptions
 from world_maker.world_maker import *
 from House import *
 
+
 def main():
 
     rectangle_house_mountain, rectangle_building, skeleton_highway, skeleton_mountain = world_maker()
 
-
     editor = Editor(buffering=True)
     buildArea = editor.getBuildArea()
 
@@ -28,19 +28,22 @@ def main():
         "garden_outline": "oak_leaves",
         "garden_floor": "grass_block"
     }
-    
+
     entranceDirection = ["N", "S", "E", "W"]
 
-    
     for houses in rectangle_house_mountain:
-        start = (houses[0][0]+buildArea.begin[0], houses[0][1], houses[0][2]+buildArea.begin[2])
+        start = (houses[0][0]+buildArea.begin[0], houses[0]
-        end = (houses[1][0]+buildArea.begin[0], houses[1][1], houses[1][2]+buildArea.begin[2])
+                 [1], houses[0][2]+buildArea.begin[2])
-        house = House(editor, start, end, entranceDirection[random.randint(0, 3)], blocks)
+        end = (houses[1][0]+buildArea.begin[0], houses[1]
+               [1], houses[1][2]+buildArea.begin[2])
+        house = House(editor, start, end,
+                      entranceDirection[random.randint(0, 3)], blocks)
         house.build()
 
+
 if __name__ == '__main__':
     main()
-    
+
 """
 from gdpc import Editor, Block, geometry, Transform
 import networks.curve as curve
@@ -74,4 +77,4 @@ geometry.placeCuboid(editor, (-5,0,-8), (25,100,25), Block("air"))
 building = Building(random_data["buildings"], [(0,0,0), (20,30,20)], baseShape, DIRECTION.EAST)
 # build it with your custom materials
 building.build(editor, ["stone_bricks","glass_pane","glass","cobblestone_wall","stone_brick_stairs","oak_planks","white_concrete","cobblestone","stone_brick_slab","iron_bars"])
-"""
+"""

networks/curve.py

@@ -1,42 +0,0 @@
-import numpy as np
-from scipy import interpolate
-
-
-class Curve:
-    def __init__(self, target_points):
-        # list of points to [(x1, y1, z1), (...), ...]
-        self.target_points = target_points
-        self.computed_points = []
-
-    def compute_curve(self, resolution=40):
-        """
-        Fill self.computed_points with a list of points that approximate a smooth curve following self.target_points.
-
-        https://stackoverflow.com/questions/18962175/spline-interpolation-coefficients-of-a-line-curve-in-3d-space
-
-        Args:
-            points (np.array): Points where the curve should pass in order.
-            resolution (int, optional): Total number of points to compute. Defaults to 40.
-        """
-        # Remove duplicates. Curve can't intersect itself
-        points = tuple(map(tuple, np.array(self.target_points)))
-        points = sorted(set(points), key=points.index)
-
-        # Change coordinates structure to (x1, x2, x3, ...), (y1, y2, y3, ...) (z1, z2, z3, ...)
-        coords = np.array(points, dtype=np.float32)
-        x = coords[:, 0]
-        y = coords[:, 1]
-        z = coords[:, 2]
-
-        # Compute
-        tck, u = interpolate.splprep([x, y, z], s=2, k=2)
-        x_knots, y_knots, z_knots = interpolate.splev(tck[0], tck)
-        u_fine = np.linspace(0, 1, resolution)
-        x_fine, y_fine, z_fine = interpolate.splev(u_fine, tck)
-
-        x_rounded = np.round(x_fine).astype(int)
-        y_rounded = np.round(y_fine).astype(int)
-        z_rounded = np.round(z_fine).astype(int)
-
-        self.computed_points = [(x, y, z) for x, y, z in zip(
-            x_rounded, y_rounded, z_rounded)]

networks/geometry/Circle.py

@@ -0,0 +1,142 @@
+from math import cos, pi, sin
+from typing import List
+
+import numpy as np
+
+from networks.geometry.Point2D import Point2D
+
+
+class Circle:
+    def __init__(self, center: Point2D):
+        self.center = center
+
+        self.radius = None
+        self.points: List[Point2D] = []
+
+        self.inner = None
+        self.outer = None
+        self.points_thick: List[Point2D] = []
+
+        self.spaced_radius = None
+        self.spaced_points: List[Point2D] = []
+
+    def __repr__(self):
+        return f"Circle(center: {self.center}, radius: {self.radius}, spaced_radius: {self.spaced_radius}, inner: {self.inner}, outer: {self.outer})"
+
+    def circle(self, radius: int) -> List[Point2D]:
+        self.radius = radius
+        center = self.center.copy()
+
+        x = -radius
+        y = 0
+        error = 2-2*radius
+        while (True):
+            self.points.append(Point2D(center.x-x, center.y+y))
+            self.points.append(Point2D(center.x-y, center.y-x))
+            self.points.append(Point2D(center.x+x, center.y-y))
+            self.points.append(Point2D(center.x+y, center.y+x))
+            r = error
+            if (r <= y):
+                y += 1
+                error += y*2+1
+            if (r > x or error > y):
+                x += 1
+                error += x*2+1
+            if (x < 0):
+                continue
+            else:
+                break
+        return self.points
+
+    def circle_thick(self, inner: int, outer: int) -> List[Point2D]:
+        """Compute discrete value of a 2d-circle with thickness. 
+
+        From: https://stackoverflow.com/questions/27755514/circle-with-thickness-drawing-algorithm
+
+        Args:
+            inner (int): The minimum radius at which the disc is filled (included).
+            outer (int): The maximum radius where disc filling stops (included).
+
+        Returns:
+            list(Point2D): List of 2d-coordinates composing the surface. Note that some coordinates are redondant and are not ordered.
+
+        >>> Circle(Point2D(0, 0), 5, 10)
+        """
+
+        self.inner = inner
+        self.outer = outer
+        center = self.center.copy()
+
+        xo = outer
+        xi = inner
+
+        y = 0
+        erro = 1 - xo
+        erri = 1 - xi
+
+        while xo >= y:
+            self._x_line(center.x + xi, center.x + xo, center.y + y)
+            self._y_line(center.x + y,  center.y + xi, center.y + xo)
+            self._x_line(center.x - xo, center.x - xi, center.y + y)
+            self._y_line(center.x - y,  center.y + xi, center.y + xo)
+            self._x_line(center.x - xo, center.x - xi, center.y - y)
+            self._y_line(center.x - y,  center.y - xo, center.y - xi)
+            self._x_line(center.x + xi, center.x + xo, center.y - y)
+            self._y_line(center.x + y,  center.y - xo, center.y - xi)
+
+            y += 1
+
+            if erro < 0:
+                erro += 2 * y + 1
+            else:
+                xo -= 1
+                erro += 2 * (y - xo + 1)
+
+            if y > inner:
+                xi = y
+            else:
+                if erri < 0:
+                    erri += 2 * y + 1
+                else:
+                    xi -= 1
+                    erri += 2 * (y - xi + 1)
+        return self.points_thick
+
+    def circle_spaced(self, number: int, radius: int) -> List[Point2D]:
+        """Get evenly spaced coordinates of the circle.
+
+        From: https://stackoverflow.com/questions/8487893/generate-all-the-points-on-the-circumference-of-a-circle
+
+        Args:
+            number (int): Number of coordinates to be returned. 
+            radius (int): Radius of the circle.
+
+        Returns:
+            list(Point2D): List of evenly spaced 2d-coordinates forming the circle.
+        """
+        self.spaced_radius = radius
+        center = self.center
+
+        self.spaced_points = [
+            Point2D(round(cos(2 * pi / number * i) * radius),
+                    round(sin(2 * pi / number * i) * radius))
+            for i in range(0, number + 1)
+        ]
+
+        for i in range(len(self.spaced_points)):
+            current_point = Point2D(
+                self.spaced_points[i].x + center.x,
+                self.spaced_points[i].y + center.y
+            ).round()
+            self.spaced_points[i] = current_point
+        return self.spaced_points
+
+    def _x_line(self, x1, x2, y):
+        while x1 <= x2:
+            self.points_thick.append(Point2D(x1, y))
+            x1 += 1
+
+    def _y_line(self, x, y1, y2):
+        while y1 <= y2:
+            self.points_thick.append(Point2D(x, y1))
+            y1 += 1

networks/geometry/Point2D.py

@@ -0,0 +1,255 @@
+from math import atan2, sqrt
+from typing import List, Union
+
+import numpy as np
+
+from Enums import ROTATION
+
+
+class Point2D:
+    def __init__(self, x: int, y: int):
+        self.x = x
+        self.y = y
+        self.coordinates = (self.x, self.y)
+
+    def copy(self):
+        return Point2D(self.x, self.y)
+
+    def __repr__(self):
+        return f"Point2D(x: {self.x}, y: {self.y})"
+
+    def __eq__(self, other):
+        if isinstance(other, Point2D):
+            return self.x == other.x and self.y == other.y
+        return False
+
+    def __add__(self, other):
+        if isinstance(other, np.ndarray) and other.shape == (2,):
+            return Point2D(self.x + other[0], self.y + other[1])
+        elif isinstance(other, Point2D):
+            return Point2D(self.x + other.x, self.y + other.y)
+        else:
+            raise TypeError(f"Unsupported type for addition: {type(other)}")
+
+    def __sub__(self, other):
+        if isinstance(other, np.ndarray) and other.shape == (2,):
+            return Point2D(self.x - other[0], self.y - other[1])
+        elif isinstance(other, Point2D):
+            return Point2D(self.x - other.x, self.y - other.y)
+        else:
+            raise TypeError(f"Unsupported type for subtraction: {type(other)}")
+
+    def is_in_triangle(self, xy0: "Point2D", xy1: "Point2D", xy2: "Point2D"):
+        """Returns True is the point is in a triangle defined by 3 others points.
+
+        From: https://stackoverflow.com/questions/2049582/how-to-determine-if-a-point-is-in-a-2d-triangle
+
+        Args:
+            xy0 (Type[Point2D]): Point of the triangle.
+            xy1 (Type[Point2D]): Point of the triangle.
+            xy2 (Type[Point2D]): Point of the triangle.
+
+        Returns:
+            bool: False if the point is not inside the triangle.
+
+        >>> Point2D(0, 0).is_in_triangle(Point2D(10, 10), Point2D(-10, 20), Point2D(0, -20)))
+        True
+        """
+        dx = self.x - xy0.x
+        dy = self.y - xy0.y
+
+        dx2 = xy2.x - xy0.x
+        dy2 = xy2.y - xy0.y
+        dx1 = xy1.x - xy0.x
+        dy1 = xy1.y - xy0.y
+
+        s_p = (dy2 * dx) - (dx2 * dy)
+        t_p = (dx1 * dy) - (dy1 * dx)
+        d = (dx1 * dy2) - (dy1 * dx2)
+
+        if d > 0:
+            return (s_p >= 0) and (t_p >= 0) and (s_p + t_p) <= d
+        else:
+            return (s_p <= 0) and (t_p <= 0) and (s_p + t_p) >= d
+
+    def nearest(self, points: List["Point2D"], return_index: bool = False) -> Union["Point2D", List[Union["Point2D", int]]]:
+        """Return the nearest point. If multiple nearest point, returns the first in the list.
+
+        Args:
+            points (List[Point2D]): List of the points to test.
+
+        Returns:
+            Point2D: The nearest point, and if multiple, the first in the list.
+        """
+        if return_index:
+            return min(
+                enumerate(points), key=lambda pair: self.distance(pair[1]))
+        return min(points, key=lambda point: self.distance(point))
+
+    def optimized_path(self, points: List["Point2D"]) -> List["Point2D"]:
+        """Get an optimized ordered path starting from the current point.
+
+        From: https://stackoverflow.com/questions/45829155/sort-points-in-order-to-have-a-continuous-curve-using-python
+
+        Args:
+            points (List[Point2D]): List of 2d-points. Could contain the current point.
+
+        Returns:
+            List[Point2D]: Ordered list of 2d-points starting from the current point.
+
+        >>> Point2D(-2, -5).optimized_path([Point2D(0, 0), Point2D(10, 5), Point2D(1, 3)])
+        [Point2D(x: -2, y: -5), Point2D(x: 0, y: 0), Point2D(x: 1, y: 3), Point2D(x: 10, y: 5)]
+        """
+        start = self
+        if start not in points:
+            points.append(start)
+        pass_by = points
+        path = [start]
+        pass_by.remove(start)
+        while pass_by:
+            nearest = min(pass_by, key=lambda point: point.distance(path[-1]))
+            path.append(nearest)
+            pass_by.remove(nearest)
+        return path
+
+    def sort_by_rotation(self, points: List["Point2D"], rotation: ROTATION = ROTATION.CLOCKWISE) -> List["Point2D"]:
+        """Sort points in clockwise order, starting from current point.
+
+        From: https://stackoverflow.com/questions/58377015/counterclockwise-sorting-of-x-y-data
+
+        Args:
+            points (List[Point2D]): List of 2d-points. Current point can be included here.
+            rotation (ROTATION): Can be ROTATION.CLOCKWISE or ROTATION.COUNTERCLOCKWISE. Optional. Defaults to ROTATION.CLOCKWISE.
+
+        Returns:
+            List[Point2D]: List of 2d-points.
+
+        >>> Point2D(-10, -10).sort_by_rotation([Point2D(10, 10), Point2D(-10, 10), Point2D(10, -10)])
+        [Point2D(x: -10, y: -10), Point2D(x: 10, y: -10), Point2D(x: 10, y: 10), Point2D(x: -10, y: 10)]
+        """
+        if self not in points:
+            points.append(self)
+        x, y = [], []
+        for i in range(len(points)):
+            x.append(points[i].x)
+            y.append(points[i].y)
+        x, y = np.array(x), np.array(y)
+
+        x0 = np.mean(x)
+        y0 = np.mean(y)
+
+        r = np.sqrt((x - x0) ** 2 + (y - y0) ** 2)
+
+        angles = np.where(
+            (y - y0) > 0,
+            np.arccos((x - x0) / r),
+            2 * np.pi - np.arccos((x - x0) / r),
+        )
+
+        mask = np.argsort(angles)
+
+        x_sorted = list(x[mask])
+        y_sorted = list(y[mask])
+
+        # Rearrange tuples to get the correct coordinates.
+        sorted_points = []
+        for i in range(len(points)):
+            j = 0
+            while (x_sorted[i] != points[j].x) and (y_sorted[i] != points[j].y):
+                j += 1
+            else:
+                sorted_points.append(Point2D(x_sorted[i], y_sorted[i]))
+
+        if rotation == ROTATION.CLOCKWISE:
+            sorted_points.reverse()
+            start_index = sorted_points.index(self)
+            return sorted_points[start_index:] + sorted_points[:start_index]
+        else:
+            start_index = sorted_points.index(self)
+            return sorted_points[start_index:] + sorted_points[:start_index]
+
+    def angle(self, xy1, xy2):
+        """
+        Compute angle (in degrees). Corner in current point.
+
+        From: https://stackoverflow.com/questions/13226038/calculating-angle-between-two-vectors-in-python
+
+        Args:
+            xy0 (numpy.ndarray): Points in the form of [x,y].
+            xy1 (numpy.ndarray): Points in the form of [x,y].
+            xy2 (numpy.ndarray): Points in the form of [x,y].
+
+        Returns:
+            float: Angle negative for counterclockwise angle, angle positive
+            for counterclockwise angle.
+
+        >>> Point2D(0, 0).angle(Point2D(10, 10), Point2D(0, -20))
+        -135.0
+        """
+        if xy2 is None:
+            xy2 = xy1.coordinate + np.array([1, 0])
+        v0 = np.array(xy1.coordinate) - np.array(self.coordinates)
+        v1 = np.array(xy2.coordinate) - np.array(self.coordinates)
+
+        angle = atan2(np.linalg.det([v0, v1]), np.dot(v0, v1))
+        return np.degrees(angle)
+
+    def round(self, ndigits: int = None) -> "Point2D":
+        self.x = round(self.x, ndigits)
+        self.y = round(self.y, ndigits)
+        self.coordinates = (self.x, self.y)
+        return self
+
+    def distance(self, point: "Point2D") -> float:
+        return sqrt((point.x - self.x) ** 2 + (point.y - self.y) ** 2)
+
+    # def slope(self, point: "Point2D") -> int:
+    #     try:
+    #         slope = (point.y - self.y) / (point.x - self.x)
+    #         return slope
+    #     except ZeroDivisionError:
+    #         return float('inf')
+
+    # def is_between_slopes(self, lower_slope: int, upper_slope: int) -> bool:
+    #     slope = self.slope(Point2D(0, 0))
+    #     print("sole", slope, (slope <= upper_slope), (slope >= lower_slope),
+    #           ((slope <= upper_slope) and (slope >= lower_slope)))
+    #     return ((slope <= upper_slope) and (slope >= lower_slope))
+
+    # def is_between_lines(self, point_1: "Point2D", point_2: "Point2D", point_a: "Point2D", point_b: "Point2D") -> bool:
+    #     slope_1, slope_a = point_1.slope(point_2), point_a.slope(point_b)
+    #     lower_slope, upper_slope = min(slope_1, slope_a), max(slope_1, slope_a)
+    #     print(self.is_between_slopes(lower_slope, upper_slope), "slope",
+    #           lower_slope, upper_slope, self.slope(Point2D(0, 0)))
+    #     print(self.x <= max(point_1.x, point_2.x, point_a.x, point_b.x), "x max")
+    #     print(self.x >= min(point_1.x, point_2.x, point_a.x, point_b.x), "x min")
+    #     print(self.y <= max(point_1.y, point_2.y, point_a.y, point_b.y), "y max")
+    #     print(self.y >= min(point_1.y, point_2.y, point_a.y, point_b.y), "y min")
+    #     return self.is_between_slopes(lower_slope, upper_slope) and self.x <= max(point_1.x, point_2.x, point_a.x, point_b.x) and self.x >= min(point_1.x, point_2.x, point_a.x, point_b.x) and self.y <= max(point_1.y, point_2.y, point_a.y, point_b.y) and self.y >= min(point_1.y, point_2.y, point_a.y, point_b.y)
+
+    @staticmethod
+    def collinear(p0: "Point2D", p1: "Point2D", p2: "Point2D") -> bool:
+        # https://stackoverflow.com/questions/9608148/python-script-to-determine-if-x-y-coordinates-are-colinear-getting-some-e
+        x1, y1 = p1.x - p0.x, p1.y - p0.y
+        x2, y2 = p2.x - p0.x, p2.y - p0.y
+        return abs(x1 * y2 - x2 * y1) < 1e-12
+
+    @staticmethod
+    def to_arrays(points: Union[List["Point2D"], "Point2D"]) -> Union[List[np.array], "Point2D"]:
+        if isinstance(points, list):
+            vectors = []
+            for point in points:
+                vectors.append(np.array(point.coordinates))
+            return vectors
+        else:
+            return np.array(points.coordinates)
+
+    @staticmethod
+    def from_arrays(vectors: Union[List[np.array], "Point2D"]) -> Union[List["Point2D"], "Point2D"]:
+        if isinstance(vectors, list):
+            points = []
+            for vector in vectors:
+                points.append(Point2D(vector[0], vector[1]))
+            return points
+        else:
+            return Point2D(vectors[0], vectors[1])

networks/geometry/Point3D.py

@@ -0,0 +1,124 @@
+from math import sqrt
+from typing import List, Union
+from networks.geometry.Point2D import Point2D
+
+import numpy as np
+
+
+class Point3D:
+    def __init__(self, x: int, y: int, z: int):
+        self.x = x
+        self.y = y
+        self.z = z
+        self.coordinates = (x, y, z)
+
+    def copy(self):
+        return Point3D(self.x, self.y, self.z)
+
+    def __repr__(self):
+        return f"Point3D(x: {self.x}, y: {self.y}, z: {self.z})"
+
+    def __eq__(self, other):
+        if isinstance(other, Point3D):
+            return self.x == other.x and self.y == other.y and self.z == other.z
+
+    def nearest(self, points: List["Point3D"]) -> "Point3D":
+        """Return the nearest point. If multiple nearest point, returns the first in the list.
+
+        Args:
+            points (List[Point2D]): List of the points to test.
+
+        Returns:
+            Point3D: The nearest point, and if multiple, the first in the list.
+
+        >>> Point3D(0, 0, 0).nearest((Point3D(-10, 10, 5), Point3D(10, 10, 1)))
+        Point3D(x: 10, y: 10, z: 1)
+        """
+        return min(points, key=lambda point: self.distance(point))
+
+    def optimized_path(self, points: List["Point3D"]) -> List["Point3D"]:
+        """Get an optimized ordered path starting from the current point.
+
+        From: https://stackoverflow.com/questions/45829155/sort-points-in-order-to-have-a-continuous-curve-using-python
+
+        Args:
+            points (List[Point3D]): List of 3d-points. Could contain the current point.
+
+        Returns:
+            List[Point3D]: Ordered list of 3d-points starting from the current point.
+
+        >>> Point3D(-2, -5, 6).optimized_path([Point3D(0, 0, 7), Point3D(10, 5, 1), Point3D(1, 3, 3)])
+        [Point3D(x: -2, y: -5, z: 6), Point3D(x: 0, y: 0, z: 7), Point3D(x: 1, y: 3, z: 3), Point3D(x: 10, y: 5, z: 1)]
+        """
+        start = self
+        if start not in points:
+            points.append(start)
+        pass_by = points
+        path = [start]
+        pass_by.remove(start)
+        while pass_by:
+            nearest = min(pass_by, key=lambda point: point.distance(path[-1]))
+            path.append(nearest)
+            pass_by.remove(nearest)
+        return path
+
+    def round(self, ndigits: int = None) -> "Point3D":
+        self.x = round(self.x, ndigits)
+        self.y = round(self.y, ndigits)
+        self.z = round(self.z, ndigits)
+        self.coordinates = (self.x, self.y, self.z)
+        return self
+
+    def distance(self, point: "Point3D"):
+        return sqrt((point.x - self.x) ** 2 + (point.y - self.y) ** 2 + (point.z - self.z) ** 2)
+
+    @staticmethod
+    def to_arrays(points: Union[List["Point3D"], "Point3D"]) -> Union[List[np.array], "Point3D"]:
+        if isinstance(points, list):
+            vectors = []
+            for point in points:
+                vectors.append(np.array(point.coordinates))
+            return vectors
+        else:
+            return np.array(points.coordinates)
+
+    @staticmethod
+    def from_arrays(vectors: Union[List[np.array], "Point3D"]) -> Union[List["Point3D"], "Point3D"]:
+        if isinstance(vectors, list):
+            points = []
+            for vector in vectors:
+                points.append(Point3D(vector[0], vector[1], vector[2]))
+            return points
+        else:
+            return Point3D(vectors[0], vectors[1], vectors[2])
+
+    @staticmethod
+    def to_2d(points: List["Point3D"], removed_axis: str) -> List[Point2D]:
+        points_2d = []
+        if removed_axis == 'x':
+            for i in range(len(points)):
+                points_2d.append(Point2D(points[i].y, points[i].z))
+        if removed_axis == 'y':
+            for i in range(len(points)):
+                points_2d.append(Point2D(points[i].x, points[i].z))
+        if removed_axis == 'z':
+            for i in range(len(points)):
+                points_2d.append(Point2D(points[i].x, points[i].y))
+        return points_2d
+
+    @staticmethod
+    def insert_3d(points: List[Point2D], position: str, to_insert: List[int]) -> List["Point3D"]:
+        points_3d = []
+        if position == 'x':
+            for i in range(len(points)):
+                points_3d.append(
+                    Point3D(to_insert[i], points[i].x, points[i].y))
+        if position == 'y':
+            for i in range(len(points)):
+                points_3d.append(
+                    Point3D(points[i].x, to_insert[i], points[i].y))
+        if position == 'z':
+            for i in range(len(points)):
+                points_3d.append(
+                    Point3D(points[i].x, points[i].y, to_insert[i]))
+        return points_3d

networks/geometry/Polyline.py

@@ -0,0 +1,231 @@
+from math import inf, sqrt
+from typing import List, Tuple, Union
+
+import numpy as np
+
+from networks.geometry.Circle import Circle
+from networks.geometry.Point2D import Point2D
+from networks.geometry.Segment2D import Segment2D
+
+
+class Polyline:
+    def __init__(self, points: List[Point2D]):
+        """A polyline with smooth corners, only composed of segments and circle arc.
+
+        Mathematics and algorithms behind this can be found here: https://cdr.lib.unc.edu/concern/dissertations/pz50gw814?locale=en, E2 Construction of arc roads from polylines, page 210.
+
+        Args:
+            points (List[Point2D]): List of 2d-points in order describing the polyline.
+
+        Raises:
+            ValueError: At least 4 points required.
+
+        >>> Polyline((Point2D(0, 0), Point2D(0, 10), Point2D(50, 10), Point2D(20, 20)))
+        """
+        self.output_points = points
+        self.points_array = Point2D.to_arrays(
+            self._remove_collinear_points(points))
+        self.length_polyline = len(self.points_array)
+
+        if self.length_polyline < 4:
+            raise ValueError("The list must contain at least 4 elements.")
+
+        self.vectors = [None] * self.length_polyline  # v
+        self.lengths = [0] * (self.length_polyline - 1)  # l
+        self.unit_vectors = [None] * self.length_polyline  # n
+        self.tangente = [0] * self.length_polyline  # f
+
+        # alpha, maximum radius factor
+        self.alpha_radii = [0] * self.length_polyline
+
+        # Useful outputs. In order to not break indexation, each list has the same length, even if for n points, there is n-2 radius.
+        # Lists will start and end with None.
+        self.radii = [0] * self.length_polyline  # r, list of points
+        self.centers = [None] * self.length_polyline  # c, list of points
+        # list of tuple of points (first intersection, corresponding corner, last intersection)
+        self.acrs_intersections = [None] * self.length_polyline
+        self.arcs = [[] for _ in range(self.length_polyline)]  # list of points
+        # self.bisectors = [None] * self.length_polyline
+
+        # For n points, there is n-1 segments. Last element should stays None.
+        self.segments = [None] * \
+            self.length_polyline  # list of segments
+
+        # Run procedure
+        self._compute_requirements()
+        self._compute_alpha_radii()
+
+        self._alpha_assign(0, self.length_polyline-1)
+        self.get_radii()
+        self.get_centers()
+        self.get_arcs_intersections()
+        self.get_arcs()
+        self.get_segments()
+
+        self.total_line_output = []
+        for i in range(1, self.length_polyline-1):
+            self.total_line_output.extend(self.segments[i].segment())
+            self.total_line_output.extend(self.arcs[i])
+        self.total_line_output.extend(
+            self.segments[self.length_polyline-1].segment())
+
+        self.total_line_output = self.total_line_output[0].optimized_path(
+            self.total_line_output)
+
+    def __repr__(self):
+        return str(self.alpha_radii)
+
+    def get_radii(self) -> List[Union[int]]:
+        for i in range(1, self.length_polyline-1):
+            self.radii[i] = round(self.alpha_radii[i] * self.tangente[i])
+        return self.radii
+
+    def get_centers(self) -> List[Union[Point2D, None]]:
+        for i in range(1, self.length_polyline-1):
+            bisector = (self.unit_vectors[i] - self.unit_vectors[i-1]) / (
+                np.linalg.norm(self.unit_vectors[i] - self.unit_vectors[i-1]))
+
+            array = self.points_array[i] + sqrt((self.radii[i]
+                                                 ** 2) + (self.alpha_radii[i] ** 2)) * bisector
+            self.centers[i] = Point2D(array[0], array[1]).round()
+        return self.centers
+
+    def get_arcs_intersections(self) -> List[Tuple[Point2D]]:
+        """Get arcs intersections points.
+
+        First and last elements elements of the list should be None. For n points, there are n-1 segments, and n-2 angle.
+
+        Returns:
+            list[tuple(Point2D)]: List of tuples composed - in order - of the first arc points, the corner points, the last arc points. The corresponding arc circle is inside this triangle.
+        """
+        for i in range(1, self.length_polyline-1):
+            point_1 = Point2D.from_arrays(self.points_array[i] -
+                                          self.alpha_radii[i] * self.unit_vectors[i-1])
+            point_2 = Point2D.from_arrays(self.points_array[i] +
+                                          self.alpha_radii[i] * self.unit_vectors[i])
+            self.acrs_intersections[i] = point_1.round(), Point2D.from_arrays(
+                self.points_array[i]), point_2.round()
+        return self.acrs_intersections
+
+    def get_arcs(self) -> List[Point2D]:
+        for i in range(1, self.length_polyline-1):
+            points = Circle(self.centers[i]).circle(self.radii[i])
+
+            # Better to do here than drawing circle arc inside big triangle!
+            double_point_a = Point2D.from_arrays(Point2D.to_arrays(self.acrs_intersections[i][0]) + 5 * (Point2D.to_arrays(
+                self.acrs_intersections[i][0]) - Point2D.to_arrays(self.centers[i])))
+            double_point_b = Point2D.from_arrays(Point2D.to_arrays(self.acrs_intersections[i][2]) + 5 * (Point2D.to_arrays(
+                self.acrs_intersections[i][2]) - Point2D.to_arrays(self.centers[i])))
+
+            for j in range(len(points)):
+                if points[j].is_in_triangle(double_point_a, self.centers[i], double_point_b):
+                    self.arcs[i].append(points[j])
+        return self.arcs
+
+    def get_segments(self) -> List[Segment2D]:
+        """Get the segments between the circle arcs and at the start and end.
+
+        Last list element should be None, and last usable index is -2 or self.length_polyline - 2. For n points, there are n-1 segments.
+
+        Returns:
+            list[Segment2D]: List of segments in order.
+        """
+        # Get first segment.
+        # First arc index is 1 because index 0 is None due to fix list lenght.  Is it a good choice?
+        self.segments[1] = Segment2D(Point2D.from_arrays(
+            self.points_array[0]), self.acrs_intersections[1][0])
+
+        # Get segments between arcs
+        for i in range(2, self.length_polyline - 1):
+            self.segments[i] = Segment2D(Point2D(self.acrs_intersections[i][0].x, self.acrs_intersections[i][0].y), Point2D(
+                self.acrs_intersections[i-1][-1].x, self.acrs_intersections[i-1][-1].y))
+
+        # Why -3?
+        # For n points, there are n-1 segments.
+        self.segments[-1] = Segment2D(self.acrs_intersections[-2][2], Point2D.from_arrays(
+            self.points_array[-1]))
+
+        return self.segments
+
+    def _alpha_assign(self, start_index: int, end_index: int):
+        """
+        The alpha-assign procedure assigning radii based on a polyline.
+        """
+        minimum_radius, minimum_index = inf, end_index
+
+        if start_index + 1 >= end_index:
+            return
+
+        alpha_b = min(
+            self.lengths[start_index] - self.alpha_radii[start_index], self.lengths[start_index + 1])
+        current_radius = max(self.tangente[start_index] * self.alpha_radii[start_index],
+                             self.tangente[start_index + 1] * alpha_b)  # Radius at initial segment
+
+        if current_radius < minimum_radius:
+            minimum_radius, minimum_index = current_radius, start_index
+            # 0, 8
+            alpha_low, alpha_high = self.alpha_radii[start_index], alpha_b
+
+        for i in range(start_index + 1, end_index - 1):  # Radii for internal segments
+            alpha_a, alpha_b, current_radius = self._radius_balance(i)
+
+            if current_radius < minimum_radius:
+                minimum_radius, minimum_index = current_radius, i
+                alpha_low, alpha_high = alpha_a, alpha_b
+
+        alpha_a = min(
+            self.lengths[end_index-2], self.lengths[end_index-1]-self.alpha_radii[end_index])
+
+        current_radius = max(self.tangente[end_index-1]*alpha_a, self.tangente[end_index]
+                             * self.alpha_radii[end_index])  # Radius at final segment
+
+        if current_radius < minimum_radius:
+            minimum_radius, minimum_index = current_radius, end_index - 1
+            alpha_low, alpha_high = alpha_a, self.alpha_radii[end_index]
+
+        # Assign alphas at ends of selected segment
+        self.alpha_radii[minimum_index] = alpha_low
+        self.alpha_radii[minimum_index+1] = alpha_high
+        # Recur on lower segments
+        self._alpha_assign(start_index, minimum_index)
+        # Recur on higher segments
+        self._alpha_assign(minimum_index + 1, end_index)
+
+    def _radius_balance(self, i: int):
+        """
+        Returns the radius that balances the radii on either end segement i.
+        """
+        alpha_a = min(self.lengths[i-1], (self.lengths[i] *
+                      self.tangente[i+1])/(self.tangente[i] + self.tangente[i+1]))
+        alpha_b = min(self.lengths[i+1], self.lengths[i]-alpha_a)
+        return alpha_a, alpha_b, min(self.tangente[i]*alpha_a, self.tangente[i+1]*alpha_b)
+
+    def _compute_requirements(self):
+        # Between two points, there is only one segment
+        for j in range(self.length_polyline-1):
+            self.vectors[j] = self.points_array[j+1] - self.points_array[j]
+            self.lengths[j] = np.linalg.norm(self.vectors[j])
+            self.unit_vectors[j] = self.vectors[j]/self.lengths[j]
+
+        # Between two segments, there is only one angle
+        for i in range(1, self.length_polyline-1):
+            dot = np.dot(self.unit_vectors[i], self.unit_vectors[i-1])
+            self.tangente[i] = sqrt((1+dot)/(1-dot))
+            # self.bisectors[i] = (self.unit_vectors[i]+self.unit_vectors[i-1]) / \
+            #     np.linalg.norm(self.unit_vectors[i]-self.unit_vectors[i-1])
+
+    def _compute_alpha_radii(self):
+        self.alpha_radii[0] = 0
+        self.alpha_radii[self.length_polyline-1] = 0
+
+    @staticmethod
+    def _remove_collinear_points(points):
+        output_points = [points[0]]
+
+        for i in range(1, len(points) - 1):
+            if not Point2D.collinear(
+                    points[i-1], points[i], points[i+1]):
+                output_points.append(points[i])
+
+        output_points.append(points[-1])
+        return output_points

networks/geometry/Segment2D.py

@@ -0,0 +1,247 @@
+from typing import List, Union
+
+import numpy as np
+
+from Enums import LINE_OVERLAP, LINE_THICKNESS_MODE
+from networks.geometry.Point2D import Point2D
+
+
+class Segment2D:
+    def __init__(self, start: Point2D, end: Point2D):
+        self.start = start
+        self.end = end
+        self.points: List[Point2D] = []
+        self.points_thick: List[Point2D] = []
+        self.thickness = None
+
+    def __repr__(self):
+        return str(f"Segment2D(start: {self.start}, end: {self.end}, points: {self.points})")
+
+    def segment(self, start: Point2D = None, end: Point2D = None, overlap: LINE_OVERLAP = LINE_OVERLAP.NONE, _is_computing_thickness: bool = False) -> Union[List[Point2D], None]:
+        """Modified Bresenham draw (line) with optional overlap.
+
+        From: https://github.com/ArminJo/Arduino-BlueDisplay/blob/master/src/LocalGUI/ThickLine.hpp
+
+        Args:
+            start (Point2D): Start point of the segment.
+            end (Point2D): End point of the segment.
+            overlap (LINE_OVERLAP): Overlap draws additional pixel when changing minor direction. For standard bresenham overlap, choose LINE_OVERLAP_NONE. Can also be LINE_OVERLAP_MAJOR or LINE_OVERLAP_MINOR.
+            _is_computing_thickness (bool, optionnal): Used by segment_thick. Don't touch.
+
+        >>> Segment2D(Point2D(0, 0), Point2D(10, 15))
+        """
+
+        if start is None or end is None:
+            start = self.start.copy()
+            end = self.end.copy()
+        else:
+            start = start.copy()
+            end = end.copy()
+
+        # Direction
+        delta_x = end.x - start.x
+        delta_y = end.y - start.y
+
+        if (delta_x < 0):
+            delta_x = -delta_x
+            step_x = -1
+        else:
+            step_x = +1
+
+        if (delta_y < 0):
+            delta_y = -delta_y
+            step_y = -1
+        else:
+            step_y = +1
+
+        delta_2x = 2*delta_x
+        delta_2y = 2*delta_y
+
+        self._add_points(start, _is_computing_thickness)
+
+        if (delta_x > delta_y):
+            error = delta_2y - delta_x
+            while (start.x != end.x):
+                start.x += step_x
+                if (error >= 0):
+                    if (overlap == LINE_OVERLAP.MAJOR):
+                        self._add_points(start, _is_computing_thickness)
+
+                    start.y += step_y
+                    if (overlap == LINE_OVERLAP.MINOR):
+                        self._add_points(
+                            Point2D(start.copy().x - step_x, start.copy().y), _is_computing_thickness)
+                    error -= delta_2x
+                error += delta_2y
+                self._add_points(start, _is_computing_thickness)
+        else:
+            error = delta_2x - delta_y
+            while (start.y != end.y):
+                start.y += step_y
+                if (error >= 0):
+                    if (overlap == LINE_OVERLAP.MAJOR):
+                        self._add_points(start, _is_computing_thickness)
+
+                    start.x += step_x
+                    if (overlap == LINE_OVERLAP.MINOR):
+                        self._add_points(
+                            Point2D(start.copy().x, start.copy().y - step_y), _is_computing_thickness)
+                    error -= delta_2y
+                error += delta_2x
+                self._add_points(start, _is_computing_thickness)
+
+        if not _is_computing_thickness:
+            return self.points
+        return None
+
+    def segment_thick(self, thickness: int, thickness_mode: LINE_THICKNESS_MODE) -> List[Point2D]:
+        """Bresenham with thickness.
+
+        From: https://github.com/ArminJo/Arduino-BlueDisplay/blob/master/src/LocalGUI/ThickLine.hpp
+        Murphy's Modified Bresenham algorithm : http://zoo.co.uk/murphy/thickline/index.html
+
+        Args:
+            start (Point2D): Start point of the segment.
+            end (Point2D): End point of the segment.
+            thickness (int): Total width of the surface. Placement relative to the original segment depends on thickness_mode.
+            thickness_mode (LINE_THICKNESS_MODE): Can be one of LINE_THICKNESS_MIDDLE, LINE_THICKNESS_DRAW_CLOCKWISE, LINE_THICKNESS_DRAW_COUNTERCLOCKWISE.
+
+        >>> self.compute_thick_segment(self.start, self.end, self.thickness, self.thickness_mode)
+        """
+        start = self.start.copy()
+        end = self.end.copy()
+
+        delta_y = end.x - start.x
+        delta_x = end.y - start.y
+
+        swap = True
+        if (delta_x < 0):
+            delta_x = -delta_x
+            step_x = -1
+            swap = not swap
+        else:
+            step_x = +1
+
+        if (delta_y < 0):
+            delta_y = -delta_y
+            step_y = -1
+            swap = not swap
+        else:
+            step_y = +1
+
+        delta_2x = 2 * delta_x
+        delta_2y = 2 * delta_y
+
+        draw_start_adjust_count = int(thickness / 2)
+        if (thickness_mode == LINE_THICKNESS_MODE.DRAW_COUNTERCLOCKWISE):
+            draw_start_adjust_count = thickness - 1
+        elif (thickness_mode == LINE_THICKNESS_MODE.DRAW_CLOCKWISE):
+            draw_start_adjust_count = 0
+
+        if (delta_x >= delta_y):
+            if swap:
+                draw_start_adjust_count = (
+                    thickness - 1) - draw_start_adjust_count
+                step_y = -step_y
+            else:
+                step_x = -step_x
+
+            error = delta_2y - delta_x
+            for i in range(draw_start_adjust_count, 0, -1):
+
+                start.x -= step_x
+                end.x -= step_x
+                if error >= 0:
+                    start.y -= step_y
+                    end.y -= step_y
+                    error -= delta_2x
+                error += delta_2x
+
+            self.segment(
+                start, end, overlap=LINE_OVERLAP.NONE, _is_computing_thickness=True)
+
+            error = delta_2x - delta_x
+            for i in range(thickness, 1, -1):
+                start.x += step_x
+                end.x += step_x
+                overlap = LINE_OVERLAP.NONE
+                if (error >= 0):
+                    start.y += step_y
+                    end.y += step_y
+                    error -= delta_2x
+                    overlap = LINE_OVERLAP.MAJOR
+                error += delta_2y
+
+                self.segment(
+                    start, end, overlap=overlap, _is_computing_thickness=True)
+
+        else:
+            if swap:
+                step_x = -step_x
+            else:
+                draw_start_adjust_count = (
+                    thickness - 1) - draw_start_adjust_count
+                step_y = -step_y
+
+            error = delta_2x - delta_y
+            for i in range(draw_start_adjust_count, 0, -1):
+                start.y -= step_y
+                end.y -= step_y
+                if (error >= 0):
+                    start.x -= step_x
+                    end.x -= step_x
+                    error -= delta_2y
+                error += delta_2x
+
+            self.segment(
+                start, end, overlap=LINE_OVERLAP.NONE, _is_computing_thickness=True)
+
+            error = delta_2x - delta_y
+            for i in range(thickness, 1, -1):
+                start.y += step_y
+                end.y += step_y
+                overlap = LINE_OVERLAP.NONE
+                if (error >= 0):
+                    start.x += step_x
+                    end.x += step_x
+                    error -= delta_2y
+                    overlap = LINE_OVERLAP.MAJOR
+                error += delta_2x
+
+                self.segment(
+                    start, end, overlap=overlap, _is_computing_thickness=True)
+
+        return self.points_thick
+
+    def perpendicular(self, distance: int) -> List[Point2D]:
+        """Compute perpendicular points from both side of the segment placed at start level.
+
+        Args:
+            distance (int): Distance bewteen the start point and the perpendicular.
+
+        Returns:
+            List[Point2D]: Two points. First one positioned on the counterclockwise side of the segment, oriented from start to end (meaning left).
+
+        >>> Segment2D(Point2D(0, 0), Point2D(10, 10)).perpendicular(10)
+        (Point2D(x: -4, y: 4), Point2D(x: 4, y: -4))
+        """
+        delta = self.start.distance(self.end)
+        dx = (self.start.x - self.end.x) / delta
+        dy = (self.start.y - self.end.y) / delta
+
+        x3 = self.start.x + (distance / 2) * dy
+        y3 = self.start.y - (distance / 2) * dx
+        x4 = self.start.x - (distance / 2) * dy
+        y4 = self.start.y + (distance / 2) * dx
+        return Point2D(x3, y3).round(), Point2D(x4, y4).round()
+
+    def middle_point(self):
+        return (np.round((self.start.x + self.end.x) / 2.0).astype(int),
+                np.round((self.start.y + self.end.y) / 2.0).astype(int),
+                )
+
+    def _add_points(self, points, is_computing_thickness):
+        if is_computing_thickness:
+            self.points_thick.append(points.copy())
+        else:
+            self.points.append(points.copy())

networks/geometry/Segment3D.py

@@ -0,0 +1,113 @@
+from typing import List
+from Enums import LINE_OVERLAP
+from networks.geometry.Point3D import Point3D
+
+
+class Segment3D:
+    def __init__(self, start: Point3D, end: Point3D):
+        self.start = start
+        self.end = end
+        self.output_points = []
+
+    def __repr__(self):
+        return str(self.output_points)
+
+    def segment(self, overlap: bool = False):
+        """Calculate a segment between two points in 3D space. 3d Bresenham algorithm.
+
+        From: https://www.geeksforgeeks.org/bresenhams-algorithm-for-3-d-line-drawing/
+
+        Args:
+            overlap (bool, optional): If False, remove unnecessary points connecting to other points side by side, leaving only a diagonal connection. Defaults to False.
+
+        >>> Segment3D(Point3D(0, 0, 0), Point3D(10, 10, 15))
+        """
+        self.output_points.append(start.copy())
+        dx = abs(self.end.x - self.start.x)
+        dy = abs(self.end.y - self.start.y)
+        dz = abs(self.end.z - self.start.z)
+        if end.x > start.x:
+            xs = 1
+        else:
+            xs = -1
+        if end.y > start.y:
+            ys = 1
+        else:
+            ys = -1
+        if end.z > start.z:
+            zs = 1
+        else:
+            zs = -1
+
+        # Driving axis is X-axis
+        if dx >= dy and dx >= dz:
+            p1 = 2 * dy - dx
+            p2 = 2 * dz - dx
+            while start.x != end.x:
+                start.x += xs
+                self.output_points.append(start.copy())
+                if p1 >= 0:
+                    start.y += ys
+                    if not overlap:
+                        if self.output_points[-1].y != start.y:
+                            self.output_points.append(start.copy())
+                    p1 -= 2 * dx
+                if p2 >= 0:
+                    start.z += zs
+                    if not overlap:
+                        if self.output_points[-1].z != start.z:
+                            self.output_points.append(start.copy())
+                    p2 -= 2 * dx
+                p1 += 2 * dy
+                p2 += 2 * dz
+
+        # Driving axis is Y-axis
+        elif dy >= dx and dy >= dz:
+            p1 = 2 * dx - dy
+            p2 = 2 * dz - dy
+            while start.y != end.y:
+                start.y += ys
+                self.output_points.append(start.copy())
+                if p1 >= 0:
+                    start.x += xs
+                    if not overlap:
+                        if self.output_points[-1].x != start.x:
+                            self.output_points.append(start.copy())
+                    p1 -= 2 * dy
+                if p2 >= 0:
+                    start.z += zs
+                    if not overlap:
+                        if self.output_points[-1].z != start.z:
+                            self.output_points.append(start.copy())
+                    p2 -= 2 * dy
+                p1 += 2 * dx
+                p2 += 2 * dz
+
+        # Driving axis is Z-axis
+        else:
+            p1 = 2 * dy - dz
+            p2 = 2 * dx - dz
+            while start.z != end.z:
+                start.z += zs
+                self.output_points.append(start.copy())
+                if p1 >= 0:
+                    start.y += ys
+                    if not overlap:
+                        if self.output_points[-1].y != start.y:
+                            self.output_points.append(start.copy())
+                    p1 -= 2 * dz
+                if p2 >= 0:
+                    start.x += xs
+                    if not overlap:
+                        if self.output_points[-1].x != start.x:
+                            self.output_points.append(start.copy())
+                    p2 -= 2 * dz
+                p1 += 2 * dy
+                p2 += 2 * dx
+        return self.output_points
+
+    def middle_point(self):
+        return (np.round((self.start.x + self.end.x) / 2.0).astype(int),
+                np.round((self.start.y + self.end.y) / 2.0).astype(int),
+                np.round((self.start.z + self.end.z) / 2.0).astype(int),
+                )

networks/geometry/Strip.py

@@ -0,0 +1,58 @@
+import networks.geometry.curve_tools as curve_tools
+import networks.geometry.segment_tools as segment_tools
+import numpy as np
+
+
+class Strip:
+    def __init__(self, points, reshape=True, spacing_distance=10):
+        self.points = np.array(points)
+        if reshape:
+            self.resolution, self.length = curve_tools.resolution_distance(
+                self.points, spacing_distance=spacing_distance)
+            self.curve = curve_tools.curve(self.points, self.resolution)
+        else:  # Point can also be given already in curved form
+            self.curve = self.points
+
+    def compute_curvature(self):
+        self.curvature = curve_tools.curvature(self.curve)
+
+    def compute_surface_perpendicular(self, width, normals):
+        self.offset_left = curve_tools.offset(self.curve, width/2, normals)
+        self.offset_right = curve_tools.offset(self.curve, -width/2, normals)
+        self.perpendicular_segment = []
+
+        for i in range(len(self.offset_left)):
+            self.perpendicular_segment.append(segment_tools.discrete_segment(
+                self.offset_left[i], self.offset_right[i], pixel_perfect=False))
+
+        self.surface = []
+
+        for i in range(len(self.perpendicular_segment)-1):
+            self.surface.append([])
+            for j in range(len(self.perpendicular_segment[i])):
+                # Hypothesis
+                max_length_index = i
+                min_length_index = i+1
+                proportion = len(
+                    self.perpendicular_segment[min_length_index])/len(self.perpendicular_segment[max_length_index])
+
+                # Reverse order if wrong hypothesis
+                if proportion > 1:
+                    max_length_index = i+1
+                    min_length_index = i
+                    proportion = len(
+                        self.perpendicular_segment[min_length_index])/len(self.perpendicular_segment[max_length_index])
+
+                self.surface[i].append([])
+                for k in range(len(self.perpendicular_segment[max_length_index])-1):
+                    self.surface[i][j].append(segment_tools.discrete_segment(
+                        self.perpendicular_segment[max_length_index][k], self.perpendicular_segment[min_length_index][round(k * proportion)], pixel_perfect=False))
+
+    def compute_surface_parallel(self, inner_range, outer_range, resolution, normals):
+        self.left_side = []
+        self.right_side = []
+        for current_range in range(inner_range * resolution, outer_range * resolution):
+            self.left_side.append(curve_tools.offset(
+                self.curve, current_range/resolution, normals))
+            self.right_side.append(curve_tools.offset(
+                self.curve, -current_range/resolution, normals))

networks/geometry/curve_tools.py

@@ -0,0 +1,141 @@
+import numpy as np
+import networks.geometry.segment_tools as segment_tools
+from scipy import interpolate
+from math import sqrt
+
+
+def curve(target_points, resolution=40):
+    """
+    Returns a list of spaced points that approximate a smooth curve following target_points.
+
+    https://stackoverflow.com/questions/18962175/spline-interpolation-coefficients-of-a-line-curve-in-3d-space
+    """
+    # Remove duplicates. Curve can't intersect itself
+    points = tuple(map(tuple, np.array(target_points)))
+    points = sorted(set(points), key=points.index)
+
+    # Change coordinates structure to (x1, x2, x3, ...), (y1, y2, y3, ...) (z1, z2, z3, ...)
+    coords = np.array(points, dtype=np.float32)
+    x = coords[:, 0]
+    y = coords[:, 1]
+    z = coords[:, 2]
+
+    # Compute
+    tck, u = interpolate.splprep([x, y, z], s=3, k=2)
+    x_knots, y_knots, z_knots = interpolate.splev(tck[0], tck)
+    u_fine = np.linspace(0, 1, resolution)
+    x_fine, y_fine, z_fine = interpolate.splev(u_fine, tck)
+
+    x_rounded = np.round(x_fine).astype(int)
+    y_rounded = np.round(y_fine).astype(int)
+    z_rounded = np.round(z_fine).astype(int)
+
+    return [(x, y, z) for x, y, z in zip(
+        x_rounded, y_rounded, z_rounded)]
+
+
+def curvature(curve):
+    """Get the normal vector at each point of the given points representing the direction in wich the curve is turning.
+
+    https://stackoverflow.com/questions/28269379/curve-curvature-in-numpy
+
+    Args:
+        curve (np.array): array of points representing the curve
+
+    Returns:
+        np.array: array of points representing the normal vector at each point in curve array
+
+    >>> curvature(np.array(([0, 0, 0], [0, 0, 1], [1, 0, 1])))
+    [[ 0.92387953 0. -0.38268343]
+    [ 0.70710678 0. -0.70710678]
+    [ 0.38268343 0. -0.92387953]]
+    """
+    curve_points = np.array(curve)
+    dx_dt = np.gradient(curve_points[:, 0])
+    dy_dt = np.gradient(curve_points[:, 1])
+    dz_dt = np.gradient(curve_points[:, 2])
+    velocity = np.array([[dx_dt[i], dy_dt[i], dz_dt[i]]
+                        for i in range(dx_dt.size)])
+
+    ds_dt = np.sqrt(dx_dt * dx_dt + dy_dt * dy_dt + dz_dt * dz_dt)
+
+    tangent = np.array([1/ds_dt]).transpose() * velocity
+    tangent_x = tangent[:, 0]
+    tangent_y = tangent[:, 1]
+    tangent_z = tangent[:, 2]
+
+    deriv_tangent_x = np.gradient(tangent_x)
+    deriv_tangent_y = np.gradient(tangent_y)
+    deriv_tangent_z = np.gradient(tangent_z)
+
+    dT_dt = np.array([[deriv_tangent_x[i], deriv_tangent_y[i], deriv_tangent_z[i]]
+                     for i in range(deriv_tangent_x.size)])
+    length_dT_dt = np.sqrt(
+        deriv_tangent_x * deriv_tangent_x + deriv_tangent_y * deriv_tangent_y + deriv_tangent_z * deriv_tangent_z + 0.0001)
+
+    normal = np.array([1/length_dT_dt]).transpose() * dT_dt
+    return normal
+
+
+def offset(curve, distance, normals):
+    if len(normals) != len(curve):
+        raise ValueError(
+            'Number of normals and number of points in the curve do not match')
+
+    # Offsetting
+    offset_segments = [segment_tools.parallel(
+        (curve[i], curve[i+1]), distance, normals[i]) for i in range(len(curve) - 1)]
+
+    # Combining segments
+    combined_curve = []
+    combined_curve.append(np.round(offset_segments[0][0]).tolist())
+    for i in range(0, len(offset_segments)-1):
+        combined_curve.append(segment_tools.middle_point(
+            offset_segments[i][1], offset_segments[i+1][0]))
+    combined_curve.append(np.round(offset_segments[-1][1]).tolist())
+
+    return combined_curve
+
+
+def resolution_distance(target_points, spacing_distance):
+    length = 0
+    for i in range(len(target_points) - 1):
+        length += sqrt(
+            ((target_points[i][0] - target_points[i + 1][0]) ** 2)
+            + ((target_points[i][1] - target_points[i + 1][1]) ** 2)
+            + ((target_points[i][2] - target_points[i + 1][2]) ** 2)
+        )
+    return round(length / spacing_distance), length
+
+
+def simplify_segments(points, epsilon):
+    if len(points) < 3:
+        return points
+
+    # Find the point with the maximum distance
+    max_distance = 0
+    max_index = 0
+    end_index = len(points) - 1
+
+    for i in range(1, end_index):
+        distance = get_distance(points[i], points[0])
+        if distance > max_distance:
+            max_distance = distance
+            max_index = i
+
+    simplified_points = []
+
+    # If the maximum distance is greater than epsilon, recursively simplify
+    if max_distance > epsilon:
+        rec_results1 = simplify_segments(points[:max_index+1], epsilon)
+        rec_results2 = simplify_segments(points[max_index:], epsilon)
+
+        # Combine the simplified sub-results
+        simplified_points.extend(rec_results1[:-1])
+        simplified_points.extend(rec_results2)
+    else:
+        # The maximum distance is less than epsilon, retain the endpoints
+        simplified_points.append(points[0])
+        simplified_points.append(points[end_index])
+
+    return simplified_points

networks/geometry/point_tools.py

@@ -0,0 +1,319 @@
+from math import sqrt, cos, pi, sin
+import numpy as np
+
+
+def segments_intersection(line0, line1, full_line=True):
+    """
+    Find (or not) intersection between two lines. Works in 2d but
+    supports 3d.
+
+    https://stackoverflow.com/questions/20677795/how-do-i-compute-the-intersection-point-of-two-lines
+
+    Args:
+        line0 (tuple): Tuple of tuple of coordinates.
+        line1 (tuple): Tuple of tuple of coordinates.
+        full_line (bool, optional): True to find intersections along
+        full line - not just in the segment.
+
+    Returns:
+        tuple: Coordinates (2d).
+
+    >>> segments_intersection(((0, 0), (0, 5)), ((2.5, 2.5), (-2.5, 2.5)))
+    """
+    xdiff = (line0[0][0] - line0[1][0], line1[0][0] - line1[1][0])
+    ydiff = (line0[0][-1] - line0[1][-1], line1[0][-1] - line1[1][-1])
+
+    def det(a, b):
+        return a[0] * b[-1] - a[-1] * b[0]
+
+    div = det(xdiff, ydiff)
+    if div == 0:
+        return None
+
+    d = (det(*line0), det(*line1))
+    x = det(d, xdiff) / div
+    y = det(d, ydiff) / div
+
+    if not full_line:
+        if (
+            min(line0[0][0], line0[1][0]) <= x <= max(line0[0][0], line0[1][0])
+            and min(line1[0][0], line1[1][0])
+            <= x
+            <= max(line1[0][0], line1[1][0])
+            and min(line0[0][-1], line0[1][-1])
+            <= y
+            <= max(line0[0][-1], line0[1][-1])
+            and min(line1[0][-1], line1[1][-1])
+            <= y
+            <= max(line1[0][-1], line1[1][-1])
+        ):
+            if len(line0[0]) > 2:
+                return middle_point(nearest(discrete_segment(line1[0], line1[1], pixel_perfect=True), (x, y)), nearest(discrete_segment(line0[0], line0[1], pixel_perfect=True), (x, y)))
+            else:
+                return x, y
+        else:
+            return None
+    else:
+        if len(line0[0]) > 2:
+            return middle_point(nearest(discrete_segment(line1[0], line1[1], pixel_perfect=True), (x, y)), nearest(discrete_segment(line0[0], line0[1], pixel_perfect=True), (x, y)))
+        else:
+            return x, y
+
+
+def circle_segment_intersection(
+    circle_center, circle_radius, xy0, xy1, full_line=True, tangent_tolerance=1e-9
+):
+    """
+    Find the points at which a circle intersects a line-segment. This
+    can happen at 0, 1, or 2 points. Works in 2d but supports 3d.
+
+    https://stackoverflow.com/questions/30844482/what-is-most-efficient-way-to-find-the-intersection-of-a-line-and-a-circle-in-py
+    Note: We follow: http://mathworld.wolfram.com/Circle-LineIntersection.html
+
+    Args:
+        circle_center (tuple): The (x, y) location of the circle center.
+        circle_radius (int): The radius of the circle.
+        xy0 (tuple): The (x, y) location of the first point of the
+        segment.
+        xy1 ([tuple]): The (x, y) location of the second point of the
+        segment.
+        full_line (bool, optional): True to find intersections along
+        full line - not just in the segment.  False will just return
+        intersections within the segment. Defaults to True.
+        tangent_tolerance (float, optional): Numerical tolerance at which we
+        decide the intersections are close enough to consider it a
+        tangent. Defaults to 1e-9.
+
+    Returns:
+        list: A list of length 0, 1, or 2, where each element is a point
+        at which the circle intercepts a line segment (2d).
+    """
+
+    (p1x, p1y), (p2x, p2y), (cx, cy) = (
+        (xy0[0], xy0[-1]),
+        (xy1[0], xy1[-1]),
+        (circle_center[0], circle_center[-1]),
+    )
+    (x1, y1), (x2, y2) = (p1x - cx, p1y - cy), (p2x - cx, p2y - cy)
+    dx, dy = (x2 - x1), (y2 - y1)
+    dr = (dx ** 2 + dy ** 2) ** 0.5
+    big_d = x1 * y2 - x2 * y1
+    discriminant = circle_radius ** 2 * dr ** 2 - big_d ** 2
+
+    if discriminant < 0:  # No intersection between circle and line
+        return []
+    else:  # There may be 0, 1, or 2 intersections with the segment
+        intersections = [
+            (
+                cx
+                + (
+                    big_d * dy
+                    + sign * (-1 if dy < 0 else 1) * dx * discriminant ** 0.5
+                )
+                / dr ** 2,
+                cy
+                + (-big_d * dx + sign * abs(dy) * discriminant ** 0.5)
+                / dr ** 2,
+            )
+            for sign in ((1, -1) if dy < 0 else (-1, 1))
+        ]  # This makes sure the order along the segment is correct
+        if (
+            not full_line
+        ):  # If only considering the segment, filter out intersections that do not fall within the segment
+            fraction_along_segment = [
+                (xi - p1x) / dx if abs(dx) > abs(dy) else (yi - p1y) / dy
+                for xi, yi in intersections
+            ]
+            intersections = [
+                pt
+                for pt, frac in zip(intersections, fraction_along_segment)
+                if 0 <= frac <= 1
+            ]
+        if (
+            len(intersections) == 2 and abs(discriminant) <= tangent_tolerance
+        ):  # If line is tangent to circle, return just one point (as both intersections have same location)
+            return [intersections[0]]
+        else:
+            return intersections
+
+
+def curved_corner_by_distance(
+    intersection, xyz0, xyz1, distance_from_intersection, resolution, full_line=True
+):
+    # Compute the merging point on the first line
+    start_curve_point_d1 = circle_segment_intersection(
+        intersection, distance_from_intersection, xyz0, intersection, full_line
+    )[0]
+    start_curve_point_d1 = (
+        round(start_curve_point_d1[0]), nearest(discrete_segment(intersection, xyz0), (start_curve_point_d1[0], 100, start_curve_point_d1[-1]))[1], round(start_curve_point_d1[-1]))
+
+    # Compute the merging point on the second line
+    end_curve_point_d1 = circle_segment_intersection(
+        intersection, distance_from_intersection, xyz1, intersection, full_line
+    )[0]
+    end_curve_point_d1 = (
+        round(end_curve_point_d1[0]), nearest(discrete_segment(intersection, xyz1), (end_curve_point_d1[0], 100, end_curve_point_d1[-1]))[1], round(end_curve_point_d1[-1]))
+
+    # Compute the merging point on the first line
+    start_curve_point_d2 = circle_segment_intersection(
+        (intersection[0], intersection[1]), distance_from_intersection, (
+            xyz0[0], xyz0[1]), (intersection[0], intersection[1]), full_line
+    )[0]
+    start_curve_point_d2 = (
+        round(start_curve_point_d2[0]), round(start_curve_point_d2[1]), nearest(discrete_segment(intersection, xyz0), (start_curve_point_d1[0], start_curve_point_d2[-1], 100))[-1])
+
+    # Compute the merging point on the second line
+    end_curve_point_d2 = circle_segment_intersection(
+        (intersection[0], intersection[1]
+         ), distance_from_intersection, (xyz1[0], xyz1[1]), (intersection[0], intersection[1]), full_line
+    )[0]
+    end_curve_point_d2 = (
+        round(end_curve_point_d2[0]), round(end_curve_point_d2[-1]), nearest(discrete_segment(
+            intersection, xyz1), (end_curve_point_d2[0], end_curve_point_d2[-1], 100))[-1])
+
+    # Compute the intersection between perpendicular lines at the merging points
+    # Higher value for better precision
+    perpendicular0_d1 = perpendicular(
+        10e3, start_curve_point_d1, intersection)[0]
+    perpendicular0_d1 = (
+        round(perpendicular0_d1[0]), round(perpendicular0_d1[-1]))
+    perpendicular1_d1 = perpendicular(
+        10e3, end_curve_point_d1, intersection)[1]
+    perpendicular1_d1 = (
+        round(perpendicular1_d1[0]), round(perpendicular1_d1[-1]))
+
+    perpendicular0_d2 = perpendicular(
+        10e3, (start_curve_point_d1[0], start_curve_point_d1[1]), (intersection[0], intersection[1]))[0]
+    perpendicular0_d2 = (
+        round(perpendicular0_d2[0]), round(perpendicular0_d2[1]))
+    perpendicular1_d2 = perpendicular(
+        10e3, (end_curve_point_d1[0], end_curve_point_d1[1]), (intersection[0], intersection[1]))[1]
+    perpendicular1_d2 = (
+        round(perpendicular1_d2[0]), round(perpendicular1_d2[1]))
+
+    # Centers
+    center_d1 = segments_intersection(
+        (perpendicular0_d1, start_curve_point_d1), (perpendicular1_d1, end_curve_point_d1))
+    center_d1 = round(center_d1[0]), middle_point(
+        xyz0, xyz1)[1], round(center_d1[-1])
+
+    center_d2 = segments_intersection(
+        (perpendicular0_d2, (start_curve_point_d1[0], start_curve_point_d1[1])), (perpendicular1_d2, (end_curve_point_d1[0], end_curve_point_d1[1])))
+    center_d2 = round(center_d2[0]), round(center_d2[1]), middle_point(
+        xyz0, xyz1)[-1]
+
+    # Compute the curvature for indications
+    curvature_d1 = round(distance(start_curve_point_d1, center_d1))
+    curvature_d2 = round(
+        distance((start_curve_point_d1[0], start_curve_point_d1[1]), center_d2))
+
+    # Return a full discrete circle or only some points of it
+    if resolution != 0:
+        circle_data_d1 = circle_points(
+            center_d1, curvature_d1, resolution
+        )
+        circle_data_d2 = circle_points(
+            center_d2, curvature_d2, resolution
+        )
+    else:
+        circle_data_d1 = circle(center_d1, curvature_d1)[0]
+        circle_data_d2 = circle(center_d2, curvature_d2)[0]
+
+    # Find the correct points on the circle.
+    curved_corner_points_temporary_d1 = [start_curve_point_d1]
+    for point in circle_data_d1:
+        if is_in_triangle(point, intersection, start_curve_point_d1, end_curve_point_d1):
+            curved_corner_points_temporary_d1.append(point)
+    curved_corner_points_temporary_d1.append(end_curve_point_d1)
+
+    # Be sure that all the points are in correct order.
+    curve_corner_points_d1 = optimized_path(
+        curved_corner_points_temporary_d1, start_curve_point_d1)
+
+    # On the other axis
+    curved_corner_points_temporary_d2 = [
+        (start_curve_point_d1[0], start_curve_point_d1[1])]
+    for point in circle_data_d2:
+
+        if is_in_triangle(point, (intersection[0], intersection[1]), (start_curve_point_d1[0], start_curve_point_d1[1]), (end_curve_point_d1[0], end_curve_point_d1[1])):
+            curved_corner_points_temporary_d2.append(point)
+    curved_corner_points_temporary_d2.append(
+        (end_curve_point_d1[0], end_curve_point_d1[1]))
+
+    # Be sure that all the points are in correct order.
+    curve_corner_points_d2 = optimized_path(
+        curved_corner_points_temporary_d2, (start_curve_point_d1[0], start_curve_point_d1[1]))
+
+    # Determine driving axis
+    if len(curve_corner_points_d1) <= len(curve_corner_points_d2):
+        main_points = curve_corner_points_d2
+        projected_points = curve_corner_points_d1
+    else:
+        main_points = curve_corner_points_d1
+        projected_points = curve_corner_points_d2
+
+    print("Main\n")
+    print(main_points)
+    print("Projected\n")
+    print(projected_points)
+
+    curve_corner_points = []
+    for i in range(len(main_points)):
+        y = projected_points[round(
+            i * (len(projected_points)-1)/len(main_points))][-1]
+        curve_corner_points.append((round(main_points[i][0]), round(
+            y), round(main_points[i][-1])))
+    return curve_corner_points, center_d1, curvature_d1, center_d2, curvature_d2
+
+
+def curved_corner_by_curvature(
+    intersection, xyz0, xyz1, curvature_radius, resolution, full_line=True
+):
+    # 3d support limited to linear interpollation on the y axis.
+    print(xyz0, intersection, xyz1)
+    # Get the center.
+    center = segments_intersection(parallel(
+        (xyz0, intersection), -curvature_radius), parallel((xyz1, intersection), curvature_radius))
+    center = round(center[0]), round(center[-1])
+
+    # Return a full discrete circle or only some points of it.
+    if resolution != 0:
+        circle_data = circle_points(
+            center, curvature_radius, resolution
+        )
+    else:
+        circle_data = circle(center, curvature_radius)[0]
+
+    # Compute the merging point on the first line.
+    print(center, curvature_radius, xyz0, intersection)
+    start_curve_point = circle_segment_intersection(
+        center, curvature_radius, xyz0, intersection, full_line
+    )[0]
+    start_curve_point = (
+        round(start_curve_point[0]), round(start_curve_point[-1]))
+
+    # Compute the merging point on the second line.
+    end_curve_point = circle_segment_intersection(
+        center, curvature_radius, xyz1, intersection, full_line
+    )[0]
+    end_curve_point = (
+        round(end_curve_point[0]), round(end_curve_point[-1]))
+
+    # Find the correct points on the circle.
+    curved_corner_points_temporary = [start_curve_point]
+    for point in circle_data:
+        # print(point, intersection, start_curve_point, end_curve_point, is_in_triangle(
+        #     point, intersection, start_curve_point, end_curve_point))
+        # if is_in_triangle(point, intersection, start_curve_point, end_curve_point):
+        curved_corner_points_temporary.append(
+            (round(point[0]), round(point[1])))
+    curved_corner_points_temporary.append(end_curve_point)
+
+    # Be sure that all the points are in correct order.
+    curve_corner_points = optimized_path(
+        curved_corner_points_temporary, start_curve_point)
+
+    # Distance from intersection just for information
+    distance_from_intersection = round(distance(start_curve_point, center))
+    return curve_corner_points, center, distance_from_intersection, parallel(
+        (xyz0, intersection), -curvature_radius), parallel((xyz1, intersection), curvature_radius)

networks/geometry/segment_tools.py

@@ -0,0 +1,56 @@
+import numpy as np
+
+
+def parallel(segment, distance, normal=np.array([0, 1, 0])):
+    """Get parallel segment in 3D space at a distance.
+
+    Args:
+        segment (np.array, np.array): start and end points of the segement.
+        distance (int): distance between both segment. Thickness in the context of a line. Positive direction means left.
+
+    Returns:
+        (np.array(), np.array()): parallel segment.
+
+    >>> parrallel(((0, 0, 0), (0, 0, 10)), 10))
+    (array([-10.,   0.,   0.]), array([-10.,   0.,  10.]))
+    """
+    return (orthogonal(segment[0], segment[1], distance, normal), orthogonal(segment[1], segment[0], -distance, normal))
+
+
+def normalized(vector):
+    magnitude = np.linalg.norm(vector)
+    if magnitude != 0:
+        normalized_vector = vector / magnitude
+        return normalized_vector
+    else:
+        return [0, 0, 0]
+
+
+def orthogonal(origin, point, distance, normal=np.array([0, 1, 0])):
+    """Get orthogonal point from a given one at the specified distance in 3D space with normal direction.
+
+    Args:
+        origin (tuple or np.array): origin
+        point (tuple or np.array): (point-origin) makes the first vector. Only the direction is used.
+        distance (int): distance from the origin. Thickness in the context of a line. Positive direction means left.
+        normal (list or np.array, optional): second vector. Defaults to the vertical [0, 1, 0].
+
+    Raises:
+        ValueError: if vectors are not linearly independent.
+
+    Returns:
+        np.array: (x y z)
+
+    >>> orthogonal((5, 5, 5), (150, 5, 5), 10)
+    [ 5.  5. 15.]
+    """
+    vector = np.subtract(point, origin)
+    normalized_vector = normalized(vector)
+    normalized_normal = normalized(normal)
+    orthogonal = np.cross(normalized_vector, normalized_normal)
+
+    if np.array_equal(orthogonal, np.zeros((3,))):
+        raise ValueError("The input vectors are not linearly independent.")
+
+    orthogonal = np.add(np.multiply(orthogonal, distance), origin).astype(int)
+    return orthogonal

networks/roads/Road.py

@@ -0,0 +1,108 @@
+import networks.geometry.curve_tools as curve_tools
+import networks.geometry.Strip as Strip
+import networks.roads.lanes.Lane as Lane
+import networks.roads.lines.Line as Line
+import json
+import random
+
+from gdpc import Editor, Block, geometry
+
+
+class Road:
+    def __init__(self, coordinates, road_configuration):
+        self.coordinates = coordinates
+        self.road_configuration = road_configuration 	# 'road', 'highway'
+        self.width = 10  # TODO
+
+    def place_roads(self):
+        editor = Editor(buffering=True)
+
+        self.resolution, self.distance = curve_tools.resolution_distance(
+            self.coordinates, 12)
+
+        self.curve_points = curve_tools.curve(
+            self.coordinates, self.resolution)
+        self.curve_surface = Strip.Strip(self.coordinates)
+        self.curve_surface.compute_curvature()
+
+        self.curvature = []
+        for i in range(len(self.curve_surface.curvature)):
+            self.curvature.append((0, 1, 0))
+
+        with open('networks/roads/lanes/lanes.json') as lanes_materials:
+            lane_type = json.load(lanes_materials).get('classic_lane')
+
+        # for coordinate, block in surface:
+        #     editor.placeBlock(coordinate, Block(block))
+
+        with open('networks/roads/lines/lines.json') as lines_materials:
+            line_type = json.load(lines_materials).get('solid_white')
+        with open('networks/roads/lines/lines.json') as lines_materials:
+            middle_line_type = json.load(lines_materials).get('broken_white')
+
+        print(line_type, lane_type)
+
+        # for coordinate, block in surface:
+        #     editor.placeBlock(coordinate, Block(block))
+
+        lines_coordinates = []
+        middle_lines_coordinates = []
+
+        # Perpendicular
+        self.curve_surface.compute_surface_perpendicular(10, self.curvature)
+        for i in range(len(self.curve_surface.surface)):
+            for j in range(len(self.curve_surface.surface[i])):
+                for k in range(len(self.curve_surface.surface[i][j])):
+                    for l in range(len(self.curve_surface.surface[i][j][k])):
+                        editor.placeBlock(
+                            self.curve_surface.surface[i][j][k][l], Block(random.choices(
+                                list(lane_type.keys()),
+                                weights=lane_type.values(),
+                                k=1,)[0]))
+                        editor.placeBlock(
+                            (self.curve_surface.surface[i][j][k][l][0], self.curve_surface.surface[i][j][k][l][1]-1, self.curve_surface.surface[i][j][k][l][2]), Block(random.choices(
+                                list(lane_type.keys()),
+                                weights=lane_type.values(),
+                                k=1,)[0]))
+                    if k == 0 or k == len(self.curve_surface.surface[i][j])-1:
+                        lines_coordinates.extend(
+                            self.curve_surface.surface[i][j][k])
+                    if k == round((len(self.curve_surface.surface[i][j])-1)/2):
+                        middle_lines_coordinates.extend(
+                            self.curve_surface.surface[i][j][k])
+
+        line = Line.Line(lines_coordinates, line_type)
+        line.get_blocks()
+        middle_line = Line.Line(middle_lines_coordinates, middle_line_type)
+        middle_line.get_blocks()
+        for i in range(len(line.coordinates_with_blocks)):
+            editor.placeBlock(
+                line.coordinates_with_blocks[i][0], Block(line.coordinates_with_blocks[i][1]))
+        for i in range(len(middle_line.coordinates_with_blocks)):
+            editor.placeBlock(
+                middle_line.coordinates_with_blocks[i][0], Block(middle_line.coordinates_with_blocks[i][1]))
+
+# offset = curve.offset(curve_surface.curve, -9, curvature)
+# for i in range(len(offset)-1):
+#     line = segment.discrete_segment(offset[i], offset[i+1])
+#     for coordinate in line:
+#         editor.placeBlock(coordinate, Block("white_concrete"))
+
+# offset = curve.offset(curve_surface.curve, 9, curvature)
+# for i in range(len(offset)-1):
+#     line = segment.discrete_segment(offset[i], offset[i+1])
+#     for coordinate in line:
+#         editor.placeBlock(coordinate, Block("white_concrete"))
+
+# # for coordinate in curve_surface.surface:
+# #     editor.placeBlock(coordinate, Block("black_concrete"))
+
+# # for coordinate in curve_surface.curve:
+# #     editor.placeBlock(coordinate, Block("red_concrete"))
+
+# # # Parallel
+# # curve_surface.compute_surface_parallel(0, 10, 8, curvature)
+
+# # for current_range in range(len(curve_surface.left_side)):
+# #     for coordinate in curve_surface.left_side[current_range]:
+# #         editor.placeBlock(coordinate, Block("yellow_concrete"))

networks/roads/intersections/Intersection.py

@@ -0,0 +1,45 @@
+from networks.geometry.segment_tools import parallel, orthogonal
+from networks.geometry.point_tools import segments_intersection, curved_corner_by_distance, curved_corner_by_curvature
+from networks.roads import Road
+
+
+class Intersection:
+    def __init__(self, center, coordinates, Roads):
+        self.center = center
+        self.coordinates = coordinates
+        self.Roads = Roads
+        self.parallel_delimitations = []
+        self.orthogonal_delimitations = []
+        self.intersections = []
+        self.intersections_curved = []
+
+    def compute_curved_corner(self):
+        # Necessary to test nearby intersection
+        self.coordinates = sort_by_clockwise(self.coordinates)
+
+        for i, coordinate in enumerate(self.coordinates):
+            right_side, left_side = parallel((coordinate, self.center), self.Roads[i].width), parallel(
+                (coordinate, self.center), -self.Roads[i].width)
+            self.parallel_delimitations.append((right_side, left_side))
+            self.orthogonal_delimitations.append(
+                ((right_side[0], left_side[0]), (right_side[-1], left_side[-1])))
+
+        for j in range(len(self.Roads)):
+            self.intersections.append(segments_intersection(
+                self.parallel_delimitations[j][1], self.parallel_delimitations[(j+1) % len(self.Roads)][0], full_line=False))
+
+            next_parallel = tuple(self.parallel_delimitations[(
+                j+1) % len(self.Roads)][0][0]), tuple(self.parallel_delimitations[(j+1) % len(self.Roads)][0][1])
+            current_parallel = tuple(self.parallel_delimitations[j][1][0]), tuple(
+                self.parallel_delimitations[j][1][1])
+
+            intersection2d = segments_intersection(((current_parallel[0][0], current_parallel[0][-1]), (current_parallel[1][0], current_parallel[1][-1])), ((
+                next_parallel[0][0], next_parallel[0][-1]), (next_parallel[1][0], next_parallel[1][-1])), full_line=False)
+
+            intersection = (
+                round(intersection2d[0]), 75, round(intersection2d[1]))
+            self.intersections_curved.append(curved_corner_by_distance(
+                intersection, current_parallel[0], next_parallel[0], 10, 0, full_line=True))
+            # print("\n\n\nBY DISTANCE\n\n\n")
+            # self.intersections_curved.append(curved_corner_by_distance(
+            #     intersection, current_parallel[0], next_parallel[0], 10, 0, full_line=True))

networks/roads/lanes/Lane.py

@@ -0,0 +1,48 @@
+import networks.geometry.curve_tools as curve_tools
+import networks.geometry.Strip as Strip
+import networks.geometry.segment_tools as segment_tools
+import random
+
+
+class Lane:
+    def __init__(self, coordinates, lane_materials,  width):
+        self.coordinates = coordinates
+        self.width = width
+        self.lane_materials = lane_materials
+        self.surface = []
+
+    def get_surface(self):
+        resolution, distance = curve_tools.resolution_distance(
+            self.coordinates, 6)
+
+        curve_points = curve_tools.curve(self.coordinates, resolution)
+        curve_surface = Strip.Strip(self.coordinates)
+        curve_surface.compute_curvature()
+
+        # Set the road to be flat
+        normals = []
+        for i in range(len(curve_surface.curvature)):
+            normals.append((0, 1, 0))
+
+        # # Compute each line
+        # for distance in range(self.width):
+        #     offset = curve_tools.offset(curve_surface.curve, distance, normals)
+        #     for i in range(len(offset)-1):
+        #         line = segment_tools.discrete_segment(offset[i], offset[i+1])
+        #         for coordinate in line:
+        #             self.surface.append((coordinate, random.choices(
+        #                 list(self.lane_materials.keys()),
+        #                 weights=self.lane_materials.values(),
+        #                 k=1,)[0]))
+
+        curve_surface.compute_surface_perpendicular(self.width, normals)
+        for i in range(len(curve_surface.surface)):
+            for j in range(len(curve_surface.surface[i])):
+                for k in range(len(curve_surface.surface[i][j])):
+                    for l in range(len(curve_surface.surface[i][j][k])):
+                        self.surface.append((curve_surface.surface[i][j][k][l], random.choices(
+                            list(self.lane_materials.keys()),
+                            weights=self.lane_materials.values(),
+                            k=1,)[0]))
+
+        return self.surface

networks/roads/lanes/lanes.json

@@ -0,0 +1,7 @@
+{
+    "classic_lane": {"stone": 3, "andesite": 1},
+    "modern_lane": {"black_concrete": 3, "black_concrete_powder": 1},
+
+    "bike_lane_green": {"green_concrete": 3, "green_concrete_powder": 1},
+    "bike_lane_red": {"red_concrete": 3, "red_concrete_powder": 1}
+}

networks/roads/lines/Line.py

@@ -0,0 +1,35 @@
+import networks.geometry.curve_tools as curve_tools
+import networks.geometry.segment_tools as segment_tools
+import random
+
+
+class Line:
+    def __init__(self, coordinates, line_materials):
+        self.coordinates = coordinates  # Full lines coordinates, not just endpoints
+        self.line_materials = line_materials  # From lines.json
+        self.coordinates_with_blocks = []  # Output
+
+    def get_blocks(self):
+        pattern_length = 0
+        pattern_materials = []
+
+        # Create the pattern materials list with correct materials depending on the selected pattern.
+        for pattern in self.line_materials:
+            pattern_length += pattern[1]
+            for _ in range(pattern[1]):
+                pattern_materials.append(pattern[0])
+
+        pattern_iteration = 0
+        for coordinate in self.coordinates:
+            block = random.choices(
+                list(pattern_materials[pattern_iteration].keys()),
+                weights=pattern_materials[pattern_iteration].values(),
+                k=1)[0]
+            if block != 'None':
+                self.coordinates_with_blocks.append((coordinate, block))
+
+            pattern_iteration += 1
+            if pattern_iteration >= pattern_length:
+                pattern_iteration = 0
+
+        return self.coordinates_with_blocks

networks/roads/lines/lines.json

@@ -0,0 +1,11 @@
+{
+    "solid_white": [
+        [{"white_concrete": 3, "white_concrete_powder": 1}, 1]
+    ],
+    
+    "broken_white": [
+        [{"red_concrete": 1}, 3],
+        [{"green_concrete": 1}, 3],
+        [{"blue_concrete": 1}, 3]
+    ]
+}

networks/roads/road.py

@@ -1,7 +0,0 @@
-class Road:
-    def __init__(self, coordinates, road_type):
-        self.coordinates = coordinates # List of tuples (x1, y1, z1) in order
-        self.road_type = road_type 	# 'road', 'highway'
-
-    def place_roads(self):
-        pass

networks/roads/roads.json

@@ -0,0 +1,11 @@
+{
+    "high_way": [
+        {"lane": "classic_lane", "width": 5, "number": 3},
+        {"lane": "classic_divider", "width": 5, "number": 1},
+        {"lane": "classic_lane", "width": 5, "number": 3}
+    ],
+    
+    "modern_road": [
+        {"lane": "classic_lane", "width": 5, "number": 2}
+    ]
+}

networks/roads_2/Roads.py

@@ -0,0 +1,89 @@
+import json
+from typing import List
+from networks.geometry.Polyline import Polyline
+
+from networks.geometry.Point3D import Point3D
+from networks.geometry.Point2D import Point2D
+from networks.geometry.Circle import Circle
+from Enums import LINE_THICKNESS_MODE
+from gdpc import Block, Editor
+
+
+class Road:
+    def __init__(self, coordinates: List[Point3D], width: int):
+        self.coordinates = coordinates
+        self.output_block = []
+        # with open(road_configuration) as f:
+        #     self.road_configuration = json.load(f)
+        #     self.width = self.road_configuration["width"]
+        self.width = width
+        self.polyline_height = None
+
+        self.polyline = Polyline(Point3D.to_2d(coordinates, 'y'))
+        self.polyline_total_line_output = [
+            [] for _ in range(len(self.polyline.total_line_output))]
+        self.index_factor = 0
+
+        self._projection()
+        self._surface()
+
+    def _surface(self):
+        # Segments
+
+        for i in range(1, len(self.polyline.segments)):
+            if len(self.polyline.segments[i].segment()) > 2:
+                for j in range(len(self.polyline.segments[i].segment_thick(self.width, LINE_THICKNESS_MODE.MIDDLE))):
+                    # Get nearest in x,z projection
+                    nearest = self.polyline.segments[i].points_thick[j].nearest(
+                        Point3D.to_2d(self.polyline_total_line_output, removed_axis='y'), True)
+                    self.output_block.append(
+                        (Point3D.insert_3d([self.polyline.segments[i].points_thick[j]], 'y', [self.polyline_total_line_output[nearest[0]].y])[0].coordinates, Block("stone")))
+
+        for i in range(1, len(self.polyline.centers)-1):
+            # Circle
+
+            circle = Circle(self.polyline.centers[i])
+            circle.circle_thick(int(
+                (self.polyline.radii[i]-self.width/2)), int((self.polyline.radii[i]+self.width/2)-1))
+
+            # Better to do here than drawing circle arc inside big triangle!
+            double_point_a = Point2D.from_arrays(Point2D.to_arrays(self.polyline.acrs_intersections[i][0]) + 5 * (Point2D.to_arrays(
+                self.polyline.acrs_intersections[i][0]) - Point2D.to_arrays(self.polyline.centers[i])))
+            double_point_b = Point2D.from_arrays(Point2D.to_arrays(self.polyline.acrs_intersections[i][2]) + 5 * (Point2D.to_arrays(
+                self.polyline.acrs_intersections[i][2]) - Point2D.to_arrays(self.polyline.centers[i])))
+
+            for j in range(len(circle.points_thick)):
+                if circle.points_thick[j].is_in_triangle(double_point_a, self.polyline.centers[i], double_point_b):
+                    nearest = circle.points_thick[j].nearest(
+                        Point3D.to_2d(self.polyline_total_line_output, removed_axis='y'), True)
+                    self.output_block.append(
+                        (Point3D.insert_3d([circle.points_thick[j]], 'y', [
+                            self.polyline_total_line_output[nearest[0]].y])[0].coordinates, Block("white_concrete")))
+
+    def _projection(self):
+        nearest_points_to_reference = []
+        for i in range(len(self.coordinates)):
+            # nearest_points_to_reference.append(Point3D.insert_3d([Point3D.to_2d([self.coordinates[i]], 'y')[0].nearest(
+            #     self.polyline.total_line_output, return_index=True)], 'y', [self.coordinates[i].y])[0])
+            index, point = Point3D.to_2d([self.coordinates[i]], 'y')[0].nearest(
+                self.polyline.total_line_output, return_index=True)
+            nearest_points_to_reference.append(
+                Point2D(index, self.coordinates[i].y))
+
+        self.polyline_height = Polyline(nearest_points_to_reference)
+
+        self.index_factor = len(
+            self.polyline_height.total_line_output)/len(self.polyline.total_line_output)
+
+        for i in range(len(self.polyline.total_line_output)):
+            self.polyline_total_line_output[i] = Point3D(
+                self.polyline.total_line_output[i].x, self.polyline_height.total_line_output[round(i*self.index_factor)].y, self.polyline.total_line_output[i].y)
+
+        self.polyline_total_line_output = self.polyline_total_line_output[0].optimized_path(
+            self.polyline_total_line_output)
+
+    def place(self):
+        editor = Editor(buffering=True)
+        for i in range(len(self.output_block)):
+            editor.placeBlock(self.output_block[i][0],
+                              self.output_block[i][1])