diff --git a/AUVSim/__init__.py b/AUVSim/__init__.py
index 63a02eb..5c439b7 100644
--- a/AUVSim/__init__.py
+++ b/AUVSim/__init__.py
@@ -1,5 +1,5 @@
 from .auv import (
-    AUV,
-    AUVPath,
-    AUVPathGen
+    AUV as AUV,
+    AUVPath as AUVPath,
+    AUVPathGen as AUVPathGen,
 )
\ No newline at end of file
diff --git a/AUVSim/auv.py b/AUVSim/auv.py
index 76819f7..73c3bb1 100644
--- a/AUVSim/auv.py
+++ b/AUVSim/auv.py
@@ -1,23 +1,31 @@
-from Raytrace.TriangleMesh import TriangleMesh, Ray
+from Raytrace import TriangleMesh
 from Raytrace.SideScan import SideScan
 import numpy as np
 from typing import Any
+from CustomTypes import Ray, Vector3, Vector3Array
+
 class AUV:
     ideal_distance:float # Currently unused
-    start_pos:np.ndarray
-    start_facing:np.ndarray
-    def __init__(self, position:np.ndarray, facing:np.ndarray):
+    start_pos:Vector3
+    start_facing:Vector3
+    def __init__(self,
+                 position:Vector3,
+                 facing:Vector3,
+                 ideal_distance:float = -1,
+                 ):
         self.start_pos = position
         self.start_facing = facing
+        self.ideal_distance = ideal_distance
 class AUVPath:
-    positions:np.ndarray[float, Any]
-    facings:np.ndarray[float, Any]
+    positions:Vector3Array
+    facings:Vector3Array
 
-    def __init__(self, positions:np.ndarray, facings:np.ndarray):
+    def __init__(self, positions:Vector3Array, facings:Vector3Array):
         self.positions = positions
         self.facings = facings
-    def as_rays(self):
-        return [Ray(i, j) for i, j in zip(self.facings, self.positions)]
+    @property
+    def rays(self) -> list[Ray]:
+        return list(map(Ray, self.facings, self.positions))
 
 class AUVPathGen:
     mesh:TriangleMesh
@@ -50,10 +58,8 @@ def get_path(self, travel_distance:float, samples:int) -> AUVPath:
             # point in the direction that the hull is sloping
             rise = (ray_dist - prev_dist) * side_ray.direction
             if not (np.isfinite(ray_dist) and np.isfinite(prev_dist)):
-                rise = np.zeros(3)
+                rise = Vector3(0, 0, 0)
             new_facing = normalize(facings[-1] + rise)
             facings.append(new_facing)
 
-        np_positions = np.array(positions)
-        np_facings = np.array(facings)
-        return AUVPath(np_positions, np_facings)
\ No newline at end of file
+        return AUVPath(Vector3Array(positions), Vector3Array(facings))
\ No newline at end of file
diff --git a/CustomTypes/.gitignore b/CustomTypes/.gitignore
new file mode 100644
index 0000000..ba0430d
--- /dev/null
+++ b/CustomTypes/.gitignore
@@ -0,0 +1 @@
+__pycache__/
\ No newline at end of file
diff --git a/CustomTypes/__init__.py b/CustomTypes/__init__.py
new file mode 100644
index 0000000..ff9b0ea
--- /dev/null
+++ b/CustomTypes/__init__.py
@@ -0,0 +1,11 @@
+from .custom_np import CustomNP as CustomNP
+from .custom_types import (
+    Vector3 as Vector3,
+    Triangle as Triangle,
+    Ray as Ray,
+)
+from .custom_arrays import (
+    FloatArray as FloatArray,
+    TriangleArray as TriangleArray,
+    Vector3Array as Vector3Array,
+)
\ No newline at end of file
diff --git a/CustomTypes/custom_arrays.py b/CustomTypes/custom_arrays.py
new file mode 100644
index 0000000..e881f0b
--- /dev/null
+++ b/CustomTypes/custom_arrays.py
@@ -0,0 +1,54 @@
+from typing import Any, Iterable, Literal, overload
+import numpy as np
+
+from . import CustomNP, Triangle, Vector3
+
+class FloatArray(CustomNP):
+    def __new__(cls,
+                content: Iterable[float]|np.ndarray|None = None,
+                ):
+        if isinstance(content, np.ndarray):
+            if content.ndim != 1: raise ValueError(content.shape) # must be 1d
+            return content.view(cls)
+        if content is None:
+            return np.ndarray((0,), float).view(cls)
+        return np.array(content).view(cls)
+    def __getitem__(self, index):
+        if isinstance(index, slice): return super().__getitem__(index).view(type(self))
+        return super().__getitem__(index)
+
+class Vector3Array(CustomNP):
+    def __new__(cls,
+                content: Iterable[Vector3]|np.ndarray|None = None,
+                ):
+        if isinstance(content, np.ndarray):
+            if content.ndim == 1: content = content.reshape((1,) + content.shape)
+            if content.ndim != 2: raise ValueError(content.shape) # must be 1d or 2d
+            if content.shape[1]  != 3: raise ValueError(content.shape) # must be shape (n, 3)
+            return content.view(cls)
+        if content is None:
+            return np.ndarray((0, 3), float).view(cls)
+        return np.array(content).view(cls)
+    def __getitem__(self, index):
+        if isinstance(index, slice): return super().__getitem__(index).view(type(self))
+        return super().__getitem__(index)
+
+class TriangleArray(CustomNP):
+    def __new__(cls,
+                content: Iterable[Triangle]|np.ndarray|None = None,
+                ):
+        if isinstance(content, np.ndarray):
+            if content.ndim == 2: content = content.reshape((1,) + content.shape)
+            if content.ndim != 3: raise ValueError(content.shape) # must be 2d or 3d
+            if content.shape[1] != 3 or content.shape[2] != 3: raise ValueError(content.shape) # must be shape (n, 3, 3)
+            return content.view(cls)
+        if content is None:
+            return np.ndarray((0, 3, 3), float).view(cls)
+        return np.array(content).view(cls)
+    def __getitem__(self, index):
+        if isinstance(index, slice): return super().__getitem__(index).view(type(self))
+        return super().__getitem__(index)
+    @property
+    def verticies(self) -> tuple[Vector3Array, Vector3Array, Vector3Array]:
+        v0, v1, v2 = self.swapaxes(0, 1)
+        return (v0.view(Vector3Array), v1.view(Vector3Array), v2.view(Vector3Array))
\ No newline at end of file
diff --git a/CustomTypes/custom_arrays.pyi b/CustomTypes/custom_arrays.pyi
new file mode 100644
index 0000000..2b29ba9
--- /dev/null
+++ b/CustomTypes/custom_arrays.pyi
@@ -0,0 +1,53 @@
+import numpy as np
+from typing import Any, Iterable, Literal, Self, Type, TypeVar, overload
+
+from . import CustomNP, Triangle, Vector3
+
+class FloatArray(CustomNP[tuple[int, Literal[3]], np.dtype[np.float_]]):
+    @overload
+    def __new__(cls) -> FloatArray:...
+    @overload
+    def __new__(cls, content: Iterable[float]) -> FloatArray:...
+    @overload
+    def __new__(cls, content: np.ndarray[tuple[int], np.dtype[np.float_]]) -> FloatArray:...
+    
+    @overload # type: ignore
+    def __getitem__(self, index:int) -> float: ...
+    @overload
+    def __getitem__(self, index:slice) -> Self: ...
+
+class Vector3Array(CustomNP[tuple[int, Literal[3]], np.dtype[np.float_]]):
+    @overload
+    def __new__(cls) -> Vector3Array:...
+    @overload
+    def __new__(cls, content: Iterable[Vector3]) -> Vector3Array:...
+    @overload
+    def __new__(cls, content: np.ndarray[tuple[int, Literal[3]], np.dtype[np.float_]]) -> Vector3Array:...
+    
+    @overload # type: ignore
+    def __getitem__(self, index:int) -> Vector3: ...
+    @overload
+    def __getitem__(self, index:tuple[int, Literal[0, 1, 2]]) -> float: ...
+    @overload
+    def __getitem__(self, index:slice) -> Self: ...
+
+class TriangleArray(CustomNP[tuple[int, Literal[3], Literal[3]], np.dtype[np.float_]]):
+    @overload
+    def __new__(cls) -> TriangleArray:...
+    @overload
+    def __new__(cls, content: Iterable[Triangle]) -> TriangleArray:...
+    @overload
+    def __new__(cls, content: np.ndarray[tuple[int, Literal[3], Literal[3]], np.dtype[np.float_]]) -> TriangleArray:...
+    
+    @overload # type: ignore
+    def __getitem__(self, index:int) -> Triangle: ...
+    @overload
+    def __getitem__(self, index:tuple[int, Literal[0, 1, 2]]) -> Vector3: ...
+    @overload
+    def __getitem__(self, index:tuple[int, Literal[0, 1, 2], Literal[0, 1, 2]]) -> float: ...
+    @overload
+    def __getitem__(self, index:slice) -> Self: ...
+
+    @property
+    def verticies(self) -> tuple[Vector3Array, Vector3Array, Vector3Array]: ...
+
diff --git a/CustomTypes/custom_np.py b/CustomTypes/custom_np.py
new file mode 100644
index 0000000..9ab0a44
--- /dev/null
+++ b/CustomTypes/custom_np.py
@@ -0,0 +1,21 @@
+
+from typing import Generic, TypeVar, Self
+import numpy as np
+
+
+_S = TypeVar('_S', bound=tuple)
+_T = TypeVar('_T', bound=np.dtype)
+class CustomNP(np.ndarray[_S, _T]):
+    def __add__(self, other):
+        return super().__add__(other).view(type(self))
+    def __sub__(self, other):
+        return super().__sub__(other).view(type(self))
+    def __mul__(self, other):
+        return super().__mul__(other).view(type(self))
+    def __truediv__(self, other):
+        return super().__truediv__(other).view(type(self))
+    def __rtruediv__(self, other):
+        return super().__rtruediv__(other).view(type(self))
+    @property
+    def nd(self):
+        return self.view(np.ndarray)
\ No newline at end of file
diff --git a/CustomTypes/custom_np.pyi b/CustomTypes/custom_np.pyi
new file mode 100644
index 0000000..9786798
--- /dev/null
+++ b/CustomTypes/custom_np.pyi
@@ -0,0 +1,15 @@
+from typing import TypeVar, Self
+import numpy as np
+
+
+_S = TypeVar('_S', bound=tuple)
+_T = TypeVar('_T', bound=np.dtype)
+class CustomNP(np.ndarray[_S, _T]):
+    def __add__(self, other:Self) -> Self: ... # type: ignore
+    def __sub__(self, other:Self) -> Self: ... # type: ignore
+    def __mul__(self, other:Self|float) -> Self: ... # type: ignore
+    def __truediv__(self, other:Self|float) -> Self: ... # type: ignore
+    def __rtruediv__(self, other:float) -> Self: ... # type: ignore
+
+    @property
+    def nd(self) -> np.ndarray: ...
diff --git a/CustomTypes/custom_types.py b/CustomTypes/custom_types.py
new file mode 100644
index 0000000..0057c89
--- /dev/null
+++ b/CustomTypes/custom_types.py
@@ -0,0 +1,42 @@
+from typing import Any, Iterable, Literal, Self, overload
+import numpy as np
+from . import CustomNP
+
+class Vector3(CustomNP):
+    def __new__(cls,
+                x: float | tuple[float, float, float] | Self, y = 0, z = 0):
+        if isinstance(x, tuple):
+            x, y, z = x
+        elif isinstance(x, Vector3):
+            x, y, z = x
+        assert isinstance(x, (float, int, np.floating))
+        return np.array([x, y, z]).view(cls)
+    def __str__(self):
+        return f"({self[0]}, {self[1]}, {self[2]})"
+    def __repr__(self):
+        return f"Vector3({self[0]}, {self[1]}, {self[2]})"
+
+class Triangle(CustomNP):
+    def __new__(cls,
+                v0: Vector3 | tuple[Vector3, Vector3, Vector3],
+                v1: Vector3 = Vector3(0, 0, 0),
+                v2: Vector3 = Vector3(0, 0, 0),
+                ) -> 'Triangle':
+        if isinstance(v0, tuple): 
+            v0, v1, v2 = v0
+        return np.array([v0, v1, v2]).view(cls)
+    
+    def __str__(self) -> str:
+        return f"({self[0]}, {self[1]}, {self[2]})"
+    def __repr__(self) -> str:
+        return f"Vector3({self[0]}, {self[1]}, {self[2]})"
+
+class Ray:
+    direction:Vector3
+    origin:Vector3
+    def __init__(self,
+                 direction:Vector3,
+                 origin:Vector3,
+                 ) -> None:
+        self.direction = direction
+        self.origin = origin
\ No newline at end of file
diff --git a/CustomTypes/custom_types.pyi b/CustomTypes/custom_types.pyi
new file mode 100644
index 0000000..04fac9a
--- /dev/null
+++ b/CustomTypes/custom_types.pyi
@@ -0,0 +1,40 @@
+from typing import Any, Iterable, Iterator, Literal, TypeAlias, TypeVar, overload
+import numpy as np
+from . import CustomNP
+
+class Vector3(CustomNP[tuple[Literal[3]], np.dtype[np.float_]]):
+    @overload
+    def __new__(cls, x: float, y: float, z: float) -> Vector3: ...
+    @overload
+    def __new__(cls, x: tuple[float, float, float]) -> Vector3: ...
+    @overload
+    def __new__(cls, x: Vector3) -> Vector3: ...
+
+    def __getitem__(self, index:Literal[0, 1, 2]) -> float: ... # type: ignore
+
+class Triangle(CustomNP[tuple[Literal[3], Literal[3]], np.dtype[np.float_]], Iterable[Vector3]):
+    @overload
+    def __new__(cls, v0: Vector3, v1: Vector3, v2: Vector3) -> Triangle: ...
+    @overload
+    def __new__(cls, v0: tuple[Vector3, Vector3, Vector3]) -> Triangle: ...
+    @overload
+    def __new__(cls, v0: Triangle) -> Triangle: ...
+
+    @overload # type: ignore
+    def __getitem__(self, index:Literal[0, 1, 2]) -> Vector3: ...
+    @overload
+    def __getitem__(self, index:tuple[Literal[0, 1, 2], Literal[0, 1, 2]]) -> float: ...
+
+    def __add__(self, other:Triangle) -> Triangle: ... # type: ignore
+    def __sub__(self, other:Triangle) -> Triangle: ... # type: ignore
+
+    def __iter__(self) -> Iterator[Vector3]: ...
+class Ray:
+    direction:Vector3
+    origin:Vector3
+    def __init__(self,
+                 direction:Vector3,
+                 origin:Vector3,
+                 ) -> None:
+        self.direction = direction
+        self.origin = origin
\ No newline at end of file
diff --git a/Generate/.gitignore b/Generate/.gitignore
new file mode 100644
index 0000000..ba0430d
--- /dev/null
+++ b/Generate/.gitignore
@@ -0,0 +1 @@
+__pycache__/
\ No newline at end of file
diff --git a/Generate/__init__.py b/Generate/__init__.py
new file mode 100644
index 0000000..94cac17
--- /dev/null
+++ b/Generate/__init__.py
@@ -0,0 +1,13 @@
+from .auv import(
+    place_auvs_over_points as place_auvs_over_points,
+    generate_paths as generate_paths,
+    scan_paths as scan_paths,
+)
+from .hull import (
+    generate_hulls as generate_hulls,
+    generate_points_on_hull as generate_points_on_hull,
+    generate_vertical_bounds as generate_vertical_bounds,
+    load_anomalies as load_anomalies,
+    pick_anomalies as pick_anomalies,
+    place_anomalies_at_points as place_anomalies_at_points,
+)
\ No newline at end of file
diff --git a/Generate/auv.py b/Generate/auv.py
new file mode 100644
index 0000000..9a82f86
--- /dev/null
+++ b/Generate/auv.py
@@ -0,0 +1,114 @@
+from typing import Callable
+
+from CustomTypes import Vector3, Vector3Array
+from Raytrace import TriangleMesh
+from AUVSim import AUV, AUVPathGen, AUVPath
+from Raytrace.SideScan import SideScan, ScanReading, RawScanReading
+
+from .helper import get_surface_normal, normalize
+
+def place_auv_over_point(
+        point:Vector3,
+        mesh:TriangleMesh,
+        distance: float,
+        align_start_facing = True,
+        use_cw_normal:None|bool = None,
+        ) -> AUV:
+    normal = get_surface_normal(mesh, point, use_cw_normal)
+    position = point + normal * distance
+    facing = Vector3(1, 0, 0)
+    if align_start_facing:
+        facing = normalize(Vector3(-normal[1], normal[0], 0)) # tangent w/ no z
+    
+    return AUV(position, facing, distance)
+
+def place_auvs_over_points(
+        points:Vector3Array,
+        meshs:TriangleMesh|list[TriangleMesh],
+        distance: float,
+        align_start_facing = True,
+        use_cw_normal:None|bool = None,
+        ) -> list[AUV]:
+    if not isinstance(meshs, list):
+        meshs = [meshs] * len(points)
+    return [
+        place_auv_over_point(point, mesh, distance, align_start_facing, use_cw_normal)
+        for point, mesh in zip(points, meshs)
+    ]
+
+def generate_path(
+        mesh:TriangleMesh,
+        auv:AUV,
+        travel_distance:float,
+        samples:int,
+        back_steps:int|None = None,
+        ) -> AUVPath:
+    # Step 1: Backstep
+    if not (back_steps is None):
+        auv_copy = AUV(auv.start_pos, auv.start_facing)
+        auv_copy.ideal_distance = auv.ideal_distance
+        back_dist = travel_distance / (samples - 1) * back_steps
+        auv = step_back_auv(mesh, auv_copy, back_dist, back_steps)
+    # Step 2: Get path
+    return AUVPathGen(mesh, auv).get_path(travel_distance, samples)
+
+def generate_paths(
+        meshs:TriangleMesh|list[TriangleMesh],
+        auvs:list[AUV],
+        travel_distance:float,
+        samples:int,
+        back_steps:int|None|list[int|None] = None,
+        ) -> list[AUVPath]:
+    if not isinstance(meshs, list):
+        meshs = [meshs] * len(auvs)
+    if not isinstance(back_steps, list): 
+        back_steps = [back_steps] * len(auvs)
+    return [
+        generate_path(mesh, auv, travel_distance, samples, bs)
+        for mesh, auv, bs in zip(meshs, auvs, back_steps)
+        ]
+
+def step_back_auv(
+        mesh:TriangleMesh,
+        auv:AUV,
+        travel_distance:float,
+        steps:int,
+        ) -> AUV:
+    auv.start_facing = auv.start_facing * -1
+    path_gen = AUVPathGen(mesh, auv)
+    back_start = path_gen.get_path(travel_distance, steps + 1).rays[-1]
+    auv.start_facing = auv.start_facing * -1 # undo changes
+
+    return AUV(back_start.origin, back_start.direction * -1, auv.ideal_distance)
+
+def scan_path(
+        scanner:SideScan,
+        path: AUVPath,
+        min_angle:float,
+        max_angle:float,
+        angle_reselution:int,
+        silent:bool = False,
+        process:Callable[[RawScanReading], ScanReading] = ScanReading
+        ) -> list[ScanReading]:
+    readings:list[ScanReading] = []
+    for n, orientation in enumerate(path.rays):
+        rays = SideScan.generate_rays(orientation, min_angle, max_angle, angle_reselution)
+        if not silent:
+            print(n + 1, len(path.rays), sep='/', end='     \r')
+        readings.append(process(scanner.scan_rays(rays)))
+    return readings
+def scan_paths(
+        scanners:SideScan|list[SideScan],
+        paths: list[AUVPath],
+        min_angle:float,
+        max_angle:float,
+        angle_reselution:int,
+        silent:bool = False,
+        process:Callable[[RawScanReading], ScanReading] = ScanReading
+        ) -> list[list[ScanReading]]:
+    if not isinstance(scanners, list): 
+        scanners = [scanners] * len(paths)
+    return [
+        scan_path(scanner, path, min_angle, max_angle, angle_reselution, silent, process)
+        for scanner, path in zip(scanners, paths)
+        ]
diff --git a/Generate/helper.py b/Generate/helper.py
new file mode 100644
index 0000000..d255759
--- /dev/null
+++ b/Generate/helper.py
@@ -0,0 +1,38 @@
+import numpy as np
+
+from CustomTypes import Ray, Vector3
+from Raytrace import TriangleMesh
+
+
+def normalize(v:Vector3) -> Vector3:
+    return v / float(np.sqrt(v.dot(v)))
+
+def get_surface_normal(
+        mesh:TriangleMesh,
+        point:Vector3,
+        use_cw_normal:None|bool = None,
+        ):
+    # Step 1: Find the triangle at the point
+    # NOTE: Using the closest centroid is not necesarily correct
+    #       But its a decent aproximation in most cases
+    min_ind = find_closest_centroid_ind(mesh, point)
+    tri = mesh.triangles[min_ind]
+    if tri is None: raise ValueError(f"sample_ray could not find surface at point {point}")
+    # Step 2: Find the normal at the triangle
+    normal = np.cross(tri[1] - tri[0], tri[2] - tri[0]).view(Vector3)
+    normal = normalize(normal)
+    if use_cw_normal:
+        normal = normal * -1
+    if use_cw_normal is None: # Guess at correct normal direction
+        if (normal[1] < 0) != (point[1] < 0):
+            normal = normal * -1 # Point normal away from xz plane
+    return normal
+
+def find_closest_centroid_ind(
+        mesh:TriangleMesh,
+        point:Vector3,
+        ) -> int:
+    # NOTE: Could (probably) be sped up with scipy.spatial.KDTree
+    rel_pos = mesh.centroids - point
+    sqr_dist = np.sum(rel_pos*rel_pos, axis=-1)
+    return int(np.argmin(sqr_dist))
\ No newline at end of file
diff --git a/Generate/hull.py b/Generate/hull.py
new file mode 100644
index 0000000..55a33aa
--- /dev/null
+++ b/Generate/hull.py
@@ -0,0 +1,122 @@
+from copy import deepcopy
+from typing import overload
+
+import numpy as np
+from CustomTypes import Vector3, TriangleArray, Vector3Array
+from Raytrace import TriangleMesh, CompositeMesh, BVHMesh
+import ShipGen as SG
+import random
+
+def generate_hulls(
+        find_best:int, 
+        out_of:int,
+        x_length:float|None = None,
+        characteristics:SG.ShipCharacteristics|None = None,
+        silent:bool = False,
+        ) -> list[BVHMesh]:
+    # Step 1: Generate the hulls
+    raw_hulls = SG.generate_hulls(find_best, out_of, characteristics)
+    # Step 2: Apply scale
+    if not (x_length is None):
+        for hull in raw_hulls:
+            SG.normalize_x(hull, x_length)
+    # Step 3: Convert them to BVHs
+    triangles_list = SG.mesh_list_to_triangles(raw_hulls)
+    hulls:list[BVHMesh] = []
+    for n, triangles in enumerate(triangles_list):
+        if not silent:
+            print(f'Building BVH {n} / {len(triangles_list)}')
+        hulls.append(BVHMesh(triangles, min_node_size = 100))
+    return hulls
+
+def place_anomaly_at_point(
+        hull:TriangleMesh,
+        anomaly:TriangleMesh,
+        point:Vector3,
+        ) -> CompositeMesh:
+    anomaly = deepcopy(anomaly)
+    # TODO: A proper addition method should be implemented
+    # HACK: Will only work fully for regular TriangleMesh
+    if type(anomaly) != TriangleMesh:
+        raise NotImplementedError('place_anomaly_at_point only works for anomalies of type TriangleMesh (no subclasses)')
+    anomaly.triangles = TriangleArray(anomaly.triangles + point)
+    
+    return CompositeMesh([hull, anomaly])
+def place_anomalies_at_points(
+        hull:TriangleMesh,
+        anomalies:list[TriangleMesh],
+        points:Vector3Array,
+        ) -> list[CompositeMesh]:
+    return [
+        place_anomaly_at_point(hull, anomaly, point)
+        for anomaly, point in zip(anomalies, points)
+        ]
+
+def generate_point_on_hull_anywhere(
+        hull:TriangleMesh,
+        ) -> Vector3:
+    # NOTE: not a perfectly even distribution
+    triangle_ind = random.randint(0, len(hull.triangles))
+    triangle = hull.triangles[triangle_ind]
+    rand_x = random.random()
+    rand_y = random.random()
+    if rand_y > 1 - rand_x: # outside triangle
+        rand_x = 1 - rand_x
+        rand_y = 1 - rand_y
+    off_x = (triangle[1] - triangle[0]) * rand_x
+    off_y = (triangle[2] - triangle[0]) * rand_y
+    return triangle[0] + off_x + off_y
+def generate_point_on_hull_in(
+        hull:TriangleMesh,
+        min_bound:Vector3,
+        max_bound:Vector3,
+        attempts:int = 1000
+        ) -> Vector3:
+    for _ in range(attempts):
+        point = generate_point_on_hull_anywhere(hull)
+        if np.all(min_bound <= point) and np.all(point <= max_bound):
+            return point
+    raise RuntimeError(f'Was not able to generate valid point within {attempts} attempts')
+
+def generate_point_on_hull(
+        hull:TriangleMesh,
+        min_bound:Vector3|None = None,
+        max_bound:Vector3|None = None,
+        attempts:int = 1000
+        ) -> Vector3:
+    if min_bound is None or max_bound is None:
+        return generate_point_on_hull_anywhere(hull)
+    else:
+        return generate_point_on_hull_in(hull, min_bound, max_bound, attempts)
+def generate_points_on_hull(
+        count:int,
+        hull:TriangleMesh,
+        min_bound:Vector3|None = None,
+        max_bound:Vector3|None = None,
+        attempts:int = 1000
+        ) -> Vector3Array:
+    return Vector3Array([
+        generate_point_on_hull(hull, min_bound, max_bound, attempts) 
+        for _ in range(count)])
+    
+def generate_vertical_bounds(
+        hull:TriangleMesh,
+        cut_percentage:float,
+        ) -> tuple[Vector3, Vector3]:
+    min_z = np.min(hull.triangles.nd[:,2])
+    max_z = np.max(hull.triangles.nd[:,2])
+    height = max_z - min_z
+    min_bound = Vector3(-np.inf, -np.inf, min_z + height * cut_percentage)
+    max_bound = Vector3( np.inf,  np.inf, max_z - height * cut_percentage)
+    return min_bound, max_bound
+
+def load_anomalies(paths:str|list[str]) -> list[TriangleMesh]:
+    if isinstance(paths, str): return load_anomalies([paths])
+    return [TriangleMesh(stl_path = path) for path in paths]
+def pick_anomalies(
+        count:int,
+        anomalies:list[TriangleMesh],
+        weights:list[float],
+        ) -> list[TriangleMesh]:
+    #TODO: Implement weights
+    return [random.choice(anomalies) for _ in range(count)]
\ No newline at end of file
diff --git a/Generate/outline.ipynb b/Generate/outline.ipynb
new file mode 100644
index 0000000..a5933cd
--- /dev/null
+++ b/Generate/outline.ipynb
@@ -0,0 +1,167 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# This notebook walks through how to create items for the dataset"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from typing import Callable\n",
+    "import sys; sys.path.append('..')\n",
+    "\n",
+    "# from Generate import (\n",
+    "#     generate_hulls,\n",
+    "#     generate_paths,\n",
+    "#     place_auvs_over_points,\n",
+    "#     generate_points_on_hull,\n",
+    "#     generate_vertical_bounds,\n",
+    "#     load_anomalies,\n",
+    "#     pick_anomalies,\n",
+    "#     place_anomalies_at_points,\n",
+    "#     scan_paths,\n",
+    "#     )\n",
+    "from Generate import *\n",
+    "\n",
+    "from ShipGen import ShipCharacteristics\n",
+    "from Raytrace.SideScan import SideScan, RawScanReading, ScanReading"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "### Parameters ###\n",
+    "hull_count = 10\n",
+    "reject_hulls = 40\n",
+    "reg_readings_per_hull = 2\n",
+    "anom_readings_per_hull = 2\n",
+    "anomalies = load_anomalies() # TODO\n",
+    "ideal_distance = 1\n",
+    "scan_distance = 1\n",
+    "readings_per_scan = 100\n",
+    "rays_per_reading = 1000\n",
+    "result_reselution = 1000\n",
+    "min_angle = np.pi/4\n",
+    "max_angle = -np.pi/2\n",
+    "\n",
+    "randomize_reading:Callable[[RawScanReading], ScanReading] = ScanReading\n",
+    "\n",
+    "# Less imortant\n",
+    "x_length = 10\n",
+    "hull_characteristics:ShipCharacteristics|None = None\n",
+    "vertical_deadzone_percent = 0.1\n",
+    "smooth_dist = 0.05"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "### Generate hulls ###\n",
+    "hulls = generate_hulls(hull_count, hull_count + reject_hulls, x_length, hull_characteristics)\n",
+    "# Loop over the hulls (only showing one for the sake of example)\n",
+    "# for hull in hulls:\n",
+    "hull = hulls[0]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "### Find points on hulls ###\n",
+    "min_b, max_b = generate_vertical_bounds(hull, vertical_deadzone_percent)\n",
+    "total_points = reg_readings_per_hull + anom_readings_per_hull\n",
+    "points_combined = generate_points_on_hull(total_points, hull, min_b, max_b)\n",
+    "reg_points = points_combined[:reg_readings_per_hull]\n",
+    "anom_points = points_combined[reg_readings_per_hull:]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "### Place anomalies ###\n",
+    "anomaly_order = pick_anomalies(anom_readings_per_hull, anomalies, [])\n",
+    "anomaly_hulls = place_anomalies_at_points(hull, anomaly_order, anom_points)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "### Place auvs ###\n",
+    "reg_auvs = place_auvs_over_points(reg_points, hull, ideal_distance)\n",
+    "anom_auvs = place_auvs_over_points(anom_points, anomaly_hulls, ideal_distance)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "### Scan hulls w/wo anomalies ###\n",
+    "# TODO: add ranomization to backsteps\n",
+    "back_steps = readings_per_scan // 2\n",
+    "reg_paths = generate_paths(hull, reg_auvs, scan_distance, readings_per_scan, back_steps)\n",
+    "scanner = SideScan(hull, smooth_dist, result_reselution)\n",
+    "reg_scans = scan_paths(scanner, reg_paths, min_angle, max_angle, rays_per_reading, process = randomize_reading)\n",
+    "\n",
+    "anom_paths = generate_paths(anomaly_hulls, anom_auvs, scan_distance, readings_per_scan, back_steps)\n",
+    "anom_scanners = [SideScan(h, smooth_dist, result_reselution) for h in anomaly_hulls]\n",
+    "anom_scans = scan_paths(anom_scanners, anom_paths, min_angle, max_angle, rays_per_reading, process = randomize_reading)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "### Done ###\n",
+    "reg_scans, anom_scans"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/Raytrace/CompositeMesh.py b/Raytrace/CompositeMesh.py
deleted file mode 100644
index f1dc5b6..0000000
--- a/Raytrace/CompositeMesh.py
+++ /dev/null
@@ -1,11 +0,0 @@
-import numpy as np
-from Raytrace.TriangleMesh import TriangleMesh, Ray
-
-class CompositeMesh(TriangleMesh):
-    meshes: list[TriangleMesh]
-    def __init__(self, meshes:list[TriangleMesh]) -> None:
-        self.meshes = meshes
-        self.triangles = np.concatenate([mesh.triangles for mesh in self.meshes])
-    def raytrace(self, ray:Ray) -> float:
-        distances = [mesh.raytrace(ray) for mesh in self.meshes]
-        return min(distances, default = np.inf)
\ No newline at end of file
diff --git a/Raytrace/Examples/ExampleMultipleSonarReadings.ipynb b/Raytrace/Examples/ExampleMultipleSonarReadings.ipynb
index d5f6cb9..9d5f4e9 100644
--- a/Raytrace/Examples/ExampleMultipleSonarReadings.ipynb
+++ b/Raytrace/Examples/ExampleMultipleSonarReadings.ipynb
@@ -2,13 +2,13 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "import sys; sys.path.append('../..')\n",
-    "from Raytrace.TriangleMesh import Ray\n",
-    "from Raytrace.BVHMesh import BVHMesh\n",
+    "from CustomTypes import Ray, Vector3, Vector3Array\n",
+    "from Raytrace import BVHMesh\n",
     "from Raytrace.SideScan import SideScan, ScanReading\n",
     "import numpy as np\n",
     "import time"
@@ -34,18 +34,18 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "facing = np.array([1, 0, 0])\n",
+    "facing = Vector3(1, 0, 0)\n",
     "min_angle = 0\n",
     "max_angle = -np.pi/2\n",
     "sample_ray_count = 1000\n",
     "\n",
-    "origin_start = np.array([5, -1, 1])\n",
-    "origin_end = np.array([9, -1, 1])\n",
+    "origin_start = Vector3(5, -1, 1)\n",
+    "origin_end = Vector3(9, -1, 1)\n",
     "readings_count = 100\n",
     "# liniar interpolation between the two for this simple demo\n",
-    "origins = [origin_start * (1-i) + origin_end * i for i in np.arange(0, 1+1/(readings_count-1)/2, 1/(readings_count-1))]\n",
+    "origins = Vector3Array([origin_start * (1-i) + origin_end * i for i in np.arange(0, 1+1/(readings_count-1)/2, 1/(readings_count-1))])\n",
     "\n",
-    "orientations = [Ray(facing, origin) for origin in origins]\n",
+    "orientations = list(map(Ray, [facing] * readings_count, origins))\n",
     "\n",
     "rays_list = [SideScan.generate_rays(orientation, min_angle, max_angle, sample_ray_count) for orientation in orientations]\n",
     "\n",
diff --git a/Raytrace/Examples/ExampleSingleSonarReading.ipynb b/Raytrace/Examples/ExampleSingleSonarReading.ipynb
index 0a89788..4c4bed6 100644
--- a/Raytrace/Examples/ExampleSingleSonarReading.ipynb
+++ b/Raytrace/Examples/ExampleSingleSonarReading.ipynb
@@ -2,16 +2,17 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "import sys; sys.path.append('../..')\n",
-    "from Raytrace.TriangleMesh import Ray\n",
-    "from Raytrace.BVHMesh import BVHMesh\n",
+    "from CustomTypes import Ray\n",
+    "from Raytrace import BVHMesh\n",
     "from Raytrace.SideScan import SideScan, ScanReading\n",
     "import numpy as np\n",
-    "import time"
+    "import time\n",
+    "from CustomTypes import Vector3"
    ]
   },
   {
@@ -34,8 +35,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "origin = np.array([5, -1, 1])\n",
-    "facing = np.array([1, 0, 0])\n",
+    "origin = Vector3(5, -1, 1)\n",
+    "facing = Vector3(1, 0, 0)\n",
     "min_angle = 0\n",
     "max_angle = -np.pi/2\n",
     "sample_ray_count = 1000\n",
@@ -47,7 +48,7 @@
     "raw_reading = SideScan(mesh).scan_rays(rays)\n",
     "reading = ScanReading(raw_reading)\n",
     "\n",
-    "print('Triangles:', mesh.triangles.shape[0])\n",
+    "print('Triangles:', len(mesh.triangles))\n",
     "raw_reading.print_summary()"
    ]
   },
diff --git a/Raytrace/Meshs/__init__.py b/Raytrace/Meshs/__init__.py
new file mode 100644
index 0000000..cc039d9
--- /dev/null
+++ b/Raytrace/Meshs/__init__.py
@@ -0,0 +1,3 @@
+from .triangle_mesh import TriangleMesh as TriangleMesh
+from .bvh_mesh import BVHMesh as BVHMesh
+from .composite_mesh import CompositeMesh as CompositeMesh
\ No newline at end of file
diff --git a/Raytrace/BVHMesh.py b/Raytrace/Meshs/bvh_mesh.py
similarity index 71%
rename from Raytrace/BVHMesh.py
rename to Raytrace/Meshs/bvh_mesh.py
index 6edb52e..14dce36 100644
--- a/Raytrace/BVHMesh.py
+++ b/Raytrace/Meshs/bvh_mesh.py
@@ -1,34 +1,37 @@
-from typing import Self, Literal
+from typing import Self, Literal, Type, TypeVar
 import numpy as np
-from Raytrace.TriangleMesh import TriangleMesh, Ray
+
+from . import TriangleMesh
+from CustomTypes import Ray, TriangleArray, Vector3, Vector3Array, Triangle
 
 # Adapted from https://jacco.ompf2.com/2022/04/13/how-to-build-a-bvh-part-1-basics/
 class BVHAABB:
-    min_pos:np.ndarray
-    max_pos:np.ndarray
+    min_pos:Vector3
+    max_pos:Vector3
     def __init__(self):
-        self.min_pos = np.array([np.nan]*3)
-        self.max_pos = np.array([np.nan]*3)
+        self.min_pos = Vector3((np.nan,)*3)
+        self.max_pos = Vector3((np.nan,)*3)
     @staticmethod
-    def from_array(triangles:np.ndarray) -> 'BVHAABB':
+    def from_array(triangles:TriangleArray) -> 'BVHAABB':
         out = BVHAABB()
-        if triangles.size == 0: return out
-        triangles = triangles.reshape((-1, 3))
-        out.min_pos = triangles.min(axis = 0)
-        out.max_pos = triangles.max(axis = 0)
+        if len(triangles) == 0: return out
+        verts = triangles.reshape((-1, 3)).view(Vector3Array)
+        out.min_pos = verts.min(axis = 0).view(Vector3)
+        out.max_pos = verts.max(axis = 0).view(Vector3)
+        
         return out
     def raytrace(self, ray:Ray) -> float:
         old_error_state = np.seterr(divide='ignore')
-        t1 = (self.min_pos - ray.origin) * ray.inverse_direction
-        t2 = (self.max_pos - ray.origin) * ray.inverse_direction
+        t1 = (self.min_pos - ray.origin) / ray.direction
+        t2 = (self.max_pos - ray.origin) / ray.direction
         np.seterr(**old_error_state)
         t = np.stack([t1, t2])
         tmin:float = t.min(axis=0).max()
         tmax:float = t.max(axis=0).min()
         return tmin if tmax >= tmin and tmax > 0 else np.inf
-    def grow(self, point:np.ndarray):
-        self.min_pos = np.nanmin(np.stack([point, self.min_pos]), axis=0)
-        self.max_pos = np.nanmax(np.stack([point, self.max_pos]), axis=0)
+    def grow(self, point:Vector3):
+        self.min_pos = np.nanmin(np.array([point, self.min_pos]), axis=0).view(Vector3)
+        self.max_pos = np.nanmax(np.array([point, self.max_pos]), axis=0).view(Vector3)
     def grow_aabb(self, other:Self):
         self.grow(other.min_pos)
         self.grow(other.max_pos)
@@ -43,27 +46,25 @@ class BVHNode:
     right:Self
     start_index:int
     tri_count:int
-    def get_subarray(self, triangles:np.ndarray) -> np.ndarray:
+    _T = TypeVar('_T', TriangleArray, Vector3Array)
+    def get_subarray(self, triangles:_T) -> _T:
         return triangles[self.start_index:self.start_index + self.tri_count]
-    def update_bounds(self, triangles:np.ndarray):
+    def update_bounds(self, triangles:TriangleArray):
         self.aabb = BVHAABB.from_array(self.get_subarray(triangles))
     def is_leaf(self) -> bool:
         return self.tri_count > 0
 class BVHMesh(TriangleMesh):
-    centroids:np.ndarray
     root:BVHNode
     node_count:int
     def __init__(self, *args, min_node_size = 2, **kwargs) -> None:
         super().__init__(*args, **kwargs)
         self.build_BVH(min_node_size)
     def build_BVH(self, min_node_size = 100) -> None:
-        # calculate triangle centroids for partitioning
-        self.centroids = self.triangles.sum(axis = 1) / 3
         # assign all triangles to root node
         self.node_count = 1
         self.root = root = BVHNode()
         root.start_index = 0
-        root.tri_count = self.triangles.shape[0]
+        root.tri_count = len(self.triangles)
         if root.tri_count == 0: return
         root.update_bounds(self.triangles)
         # subdivide recursively
@@ -80,8 +81,9 @@ def subdivide(self,  node:BVHNode, min_node_size = 100) -> None:
         while i <= j:
             if self.centroids[i, axis] < split_pos: i += 1
             else:
-                self.triangles[[i,j]] = self.triangles[[j,i]]
-                self.centroids[[i,j]] = self.centroids[[j,i]]
+                self.triangles.nd[[i,j]] = self.triangles.nd[[j,i]]
+                self.centroids.nd[[i,j]] = self.centroids.nd[[j,i]]
+
                 j -= 1
         # abort split if one of the sides is empty
         left_count = i - node.start_index
@@ -103,18 +105,22 @@ def subdivide(self,  node:BVHNode, min_node_size = 100) -> None:
         self.subdivide(left)
         self.subdivide(right)
 
-    def raytrace(self, ray:Ray) -> float:
+    def raytrace_as_tri(self, ray:Ray) -> tuple[float, Triangle|None]:
         return self.BVH_raytrace(ray)
     
-    def BVH_raytrace(self, ray:Ray) -> float:
+    def BVH_raytrace(self, ray:Ray) -> tuple[float, Triangle|None]:
         node:BVHNode = self.root
         stack:list[BVHNode] = []
-        best:float = np.inf
+        best_dist:float = np.inf
+        best_tri:Triangle|None = None
         while 1:
             if node.is_leaf():
                 intersections = self.batch_triangle_ray_intersection(node.get_subarray(self.triangles), ray)
                 intersections[intersections < 0] = np.inf
-                best = min(best, np.min(intersections))
+                min_index = int(np.argmin(intersections))
+                if intersections[min_index] < best_dist:
+                    best_dist = intersections[min_index]
+                    best_tri = node.get_subarray(self.triangles)[min_index]
                 if len(stack) == 0: break
                 node = stack.pop()
                 continue
@@ -131,25 +137,27 @@ def BVH_raytrace(self, ray:Ray) -> float:
                 node = child1
                 if np.isfinite(dist2):
                     stack.append(child2)
-        return best
-    def find_best_split(self, node:BVHNode) -> tuple[float, int, float]:
+        return best_dist, best_tri
+
+    def find_best_split(self, node:BVHNode) -> tuple[float, Literal[0, 1, 2], float]:
         BINS = 8
-        axis:int = -1
+        axis:Literal[0, 1, 2] = 0
         split_pos:float = 0
         best_cost:float = np.inf
-        triangles:np.ndarray[tuple[int, Literal[3], Literal[3]], np.dtype[np.float32]] = node.get_subarray(self.triangles)
-        centroids:np.ndarray[tuple[int, Literal[3]], np.dtype[np.float32]] = node.get_subarray(self.centroids)
-        for axis_i in range(3):
-            bounds_min = float(np.min(centroids[:,axis_i]))
-            bounds_max = float(np.max(centroids[:,axis_i]))
+        triangles:TriangleArray = node.get_subarray(self.triangles)
+        centroids:Vector3Array = node.get_subarray(self.centroids)
+        for axis_i in (0, 1, 2):
+            bounds_min = float(np.min(centroids.nd[:,axis_i]))
+            bounds_max = float(np.max(centroids.nd[:,axis_i]))
             
             if bounds_min == bounds_max: continue
             # populate the bins
             scale:float = BINS / (bounds_max - bounds_min)
-            bin_idx = ((centroids[:,axis_i] - bounds_min) * scale).astype(int)
+            bin_idx = ((centroids.nd[:,axis_i] - bounds_min) * scale).astype(int)
             bin_idx[bin_idx > BINS - 1] = BINS - 1
             bin_bounds = [BVHAABB.from_array(triangles[bin_idx == n]) for n in range(BINS)]
             bin_counts = [np.count_nonzero(bin_idx == n) for n in range(BINS)]
+            
             # gather data for the 7 planes between the 8 bins
             left_area = np.zeros(BINS - 1, float)
             right_area = np.zeros(BINS - 1, float)
diff --git a/Raytrace/Meshs/composite_mesh.py b/Raytrace/Meshs/composite_mesh.py
new file mode 100644
index 0000000..6292cb3
--- /dev/null
+++ b/Raytrace/Meshs/composite_mesh.py
@@ -0,0 +1,22 @@
+from functools import cached_property
+import numpy as np
+from . import TriangleMesh
+from CustomTypes import Ray, TriangleArray, Triangle
+
+class CompositeMesh(TriangleMesh):
+    meshes: list[TriangleMesh]
+    
+    def __init__(self, meshes:list[TriangleMesh]) -> None:
+        self.meshes = meshes
+
+    def raytrace_as_tri(self, ray:Ray) -> tuple[float, Triangle|None]:
+        traces = [mesh.raytrace_as_tri(ray) for mesh in self.meshes]
+        if len(traces) == 0: return (np.inf, None)
+        return min(traces, key= lambda v: v[0])
+    
+    @property # only should be used for rendering, not calculations
+    def triangles(self) -> TriangleArray:
+        return np.concatenate([mesh.triangles for mesh in self.meshes]).view(TriangleArray)
+    @triangles.setter
+    def triangles(self, _) -> None:
+        raise AttributeError('CompositeMesh does not support writing to its triangles property')
diff --git a/Raytrace/TriangleMesh.py b/Raytrace/Meshs/triangle_mesh.py
similarity index 55%
rename from Raytrace/TriangleMesh.py
rename to Raytrace/Meshs/triangle_mesh.py
index 50be840..e8c31ab 100644
--- a/Raytrace/TriangleMesh.py
+++ b/Raytrace/Meshs/triangle_mesh.py
@@ -1,31 +1,30 @@
+from typing import Any
 import numpy as np
 import pickle
-
-class Ray:
-    direction:np.ndarray
-    origin:np.ndarray
-    inverse_direction:np.ndarray
-    def __init__(self, direction:np.ndarray, origin:np.ndarray) -> None:
-        self.direction = direction
-        self.origin = origin
-        old_error_state = np.seterr(divide='ignore', invalid='ignore')
-        self.inverse_direction = np.ones(3) / direction
-        np.seterr(**old_error_state)
+from CustomTypes import Triangle, TriangleArray, Ray, Vector3, Vector3Array, FloatArray
+import stl
 
 class TriangleMesh:
-    # An n by 3 by 3 array of floats representing a list of triangles in the form triangle_index, vertex_index, axis_index
-    triangles:np.ndarray
-
-    def __init__(self, triangles: np.ndarray|None = None, stl_path: str|None = None):
+    triangles:TriangleArray
+    def __init__(self,
+                 triangles: np.ndarray|TriangleArray|None = None,
+                 stl_path: str|None = None):
         if not stl_path is None:
-            import stl
             mesh = stl.Mesh.from_file(stl_path)
-            triangles = np.array([mesh.v0, mesh.v1, mesh.v2]).swapaxes(0, 1)
-        self.triangles = np.ndarray((0,3), float) if triangles is None else triangles
+            self.triangles = TriangleArray(np.array([mesh.v0, mesh.v1, mesh.v2]).swapaxes(0, 1))
+        elif triangles is None:
+            self.triangles = TriangleArray([])
+        elif isinstance(triangles, np.ndarray):
+            self.triangles = TriangleArray(triangles)
+        else:
+            self.triangles = triangles
     
     def raytrace(self, ray:Ray) -> float:
-        return self.array_raytrace(ray)
+        return self.raytrace_as_tri(ray)[0]
     
+    def raytrace_as_tri(self, ray:Ray) -> tuple[float, Triangle|None]:
+        return self.array_raytrace(ray)
+
     def simple_raytrace(self, ray:Ray) -> float:
         best_dist = np.inf
         
@@ -38,40 +37,41 @@ def simple_raytrace(self, ray:Ray) -> float:
             if dist < best_dist:
                 best_dist = dist
         return best_dist
-    def array_raytrace(self, ray:Ray) -> float:
+    def array_raytrace(self, ray:Ray) -> tuple[float, Triangle|None]:
         out = self.batch_triangle_ray_intersection(self.triangles, ray)
         out[out < 0] = np.inf
-
-        return np.min(out)
+        min_index = int(np.argmin(out))
+        if np.isinf(out[min_index]): return (np.inf, None) # hit nothing
+        return (out[min_index], self.triangles[int(np.argmin(out))])
 
     @staticmethod
-    def triangle_ray_intersection(triangle, ray:Ray, epsilon = 1e-10) -> float:
+    def triangle_ray_intersection(triangle:Triangle, ray:Ray, epsilon = 1e-10) -> float:
         # Translated from https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm#C++_implementation
 
         v0, v1, v2 = triangle
         edge1 = v1 - v0
         edge2 = v2 - v0
-        ray_cross_e2 = np.cross(ray.direction, edge2)
-        det = np.dot(edge1, ray_cross_e2)
+        ray_cross_e2 = np.cross(ray.direction, edge2).view(Vector3)
+        det:float = np.dot(edge1, ray_cross_e2)
 
         if -epsilon < det < epsilon:
             return -1  # This ray is parallel to this triangle.
 
         inv_det = 1.0 / det
         s = ray.origin - v0
-        u = inv_det * np.dot(s, ray_cross_e2)
+        u:float = inv_det * np.dot(s, ray_cross_e2)
 
         if u < -epsilon or 1 < u:
             return -1
 
-        s_cross_e1 = np.cross(s, edge1)
-        v = inv_det * np.dot(ray.direction, s_cross_e1)
+        s_cross_e1 = np.cross(s, edge1).view(Vector3)
+        v:float = inv_det * np.dot(ray.direction, s_cross_e1)
 
         if v < -epsilon or 1 < u + v:
             return -1
 
         # At this stage we can compute t to find out where the intersection point is on the line.
-        t = inv_det * np.dot(edge2, s_cross_e1)
+        t:float = inv_det * np.dot(edge2, s_cross_e1)
 
         if t > epsilon:  # ray intersection
             return t
@@ -79,28 +79,27 @@ def triangle_ray_intersection(triangle, ray:Ray, epsilon = 1e-10) -> float:
         else:  # This means that there is a line intersection but not a ray intersection.
             return -1
     @staticmethod
-    def batch_triangle_ray_intersection(triangle_array, ray:Ray, epsilon = 1e-10) -> np.ndarray:
-        vecdot = lambda a, b : np.sum(a*b, axis=-1)
+    def batch_triangle_ray_intersection(triangle_array:TriangleArray, ray:Ray, epsilon = 1e-10) -> np.ndarray[np.floating, Any]:
+        def vecdot(a:Vector3Array|Vector3, b:Vector3Array) -> FloatArray: return np.sum(a*b, axis=-1).view(FloatArray)
         # Translated from https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm#C++_implementation
         
-        triangle_array_comp = triangle_array.swapaxes(0, 1)
-        v0, v1, v2 = triangle_array_comp[0], triangle_array_comp[1], triangle_array_comp[2]
-        
+        v0, v1, v2 = triangle_array.verticies
+
         # np.subtract(v1, v0, out= v1)
         # np.subtract(v2, v0, out= v2)
         v1 = v1 - v0
         v2 = v2 - v0
 
-        ray_cross_e2 = np.cross(ray.direction, v2)
+        ray_cross_e2 = np.cross(ray.direction, v2).view(Vector3Array)
         det = vecdot(v1, ray_cross_e2)
         # if -epsilon < det < epsilon:
         #     return None  # This ray is parallel to this triangle.
         old_error_state = np.seterr(divide='ignore', invalid='ignore')
         inv_det = 1.0 / det
-        s = ray.origin - v0
+        s = (ray.origin - v0).view(Vector3Array)
         u = inv_det * vecdot(s, ray_cross_e2)
 
-        s_cross_e1 = np.cross(s, v1)
+        s_cross_e1 = np.cross(s, v1).view(Vector3Array)
         v = inv_det * vecdot(ray.direction, s_cross_e1)
 
         valid = u < -epsilon
@@ -109,7 +108,7 @@ def batch_triangle_ray_intersection(triangle_array, ray:Ray, epsilon = 1e-10) ->
         np.logical_or(valid, 1 < u + v, out=valid)
         np.invert(valid, out=valid)
 
-        t = np.empty(triangle_array.shape[0], dtype=np.float64)
+        t = np.empty(len(triangle_array), dtype=float)
         t.fill(-1)
 
         np.multiply(inv_det, vecdot(v2, s_cross_e1), out = t, where = valid)
@@ -117,6 +116,16 @@ def batch_triangle_ray_intersection(triangle_array, ray:Ray, epsilon = 1e-10) ->
         np.seterr(**old_error_state)
         return t
 
+    _centroids:Vector3Array|None = None
+    @property
+    def centroids(self) -> Vector3Array:
+        if self._centroids is None:
+            self._centroids = self.triangles.sum(axis = 1).view(Vector3Array) / 3
+        return self._centroids
+    @centroids.setter
+    def centroids(self, value:Vector3Array) -> None:
+        self._centroids = value
+    
     def save(self, filename, override = False) -> None:
         with open(filename, 'wb' if override else 'xb') as file:
             pickle.dump(self, file)
@@ -126,4 +135,4 @@ def load(filename) -> 'TriangleMesh':
             obj = pickle.load(file)
             if isinstance(obj, TriangleMesh):
                 return obj
-            raise TypeError(f'The object saved in {filename} is type {type(obj)} not TriangleMesh')
\ No newline at end of file
+            raise TypeError(f'The object saved in {filename} is type {type(obj)} not TriangleMesh')
diff --git a/Raytrace/PlotRays.py b/Raytrace/PlotRays.py
index 6fa0dd2..a28dccc 100644
--- a/Raytrace/PlotRays.py
+++ b/Raytrace/PlotRays.py
@@ -1,5 +1,5 @@
-from Raytrace.TriangleMesh import TriangleMesh
-from Raytrace.SideScan import ScanReading, RawScanReading
+from . import TriangleMesh
+from .SideScan import ScanReading, RawScanReading
 import numpy as np
 try:
     import plotly.graph_objs as go # type: ignore
diff --git a/Raytrace/SideScan.py b/Raytrace/SideScan.py
index ac14c99..1900053 100644
--- a/Raytrace/SideScan.py
+++ b/Raytrace/SideScan.py
@@ -1,4 +1,5 @@
-from Raytrace.TriangleMesh import TriangleMesh, Ray
+from . import TriangleMesh
+from CustomTypes import Ray, Vector3, Vector3Array
 from typing import Self
 import numpy as np
 import time
@@ -34,7 +35,7 @@ def combine_min(self, other:Self) -> Self:
         new_distances[other_smaller] = other.distances[other_smaller]
         
         # HACK: new rays should be generated based on which distance was closer
-        new_rays = [Ray(i, j) for i, j in zip(self.directions, self.origins)] 
+        new_rays = list(map(Ray, self.directions, self.origins))
         return type(self)(new_distances, new_rays)
 
     def print_summary(self) -> None:
@@ -130,13 +131,14 @@ def generate_rays(orientation:Ray, min_angle:float, max_angle:float, angle_resel
         angle_step = (max_angle-min_angle)/(angle_reselution-1)
         angles = np.arange(min_angle, max_angle + angle_step/2, angle_step)
         origin = orientation.origin
-        pitch = np.arctan2(orientation.direction[2],np.sqrt(np.dot(orientation.direction[0:2], orientation.direction[0:2])))
-        yaw = np.arctan2(orientation.direction[1], orientation.direction[0]) # yaw = 0 => facing due +x
+        pitch:float = np.arctan2(orientation.direction[2],np.sqrt(np.dot(orientation.direction.nd[0:2], orientation.direction.nd[0:2])))
+        yaw:float = np.arctan2(orientation.direction[1], orientation.direction[0]) # yaw = 0 => facing due +x
         
         # Precomputed rotation matrix
         sin_comp = np.array([-np.sin(pitch) * np.cos(yaw), np.sin(pitch) * np.sin(yaw), np.cos(pitch)])
         cos_comp = np.array([np.sin(yaw), np.cos(yaw), 0])
 
-        return [Ray(np.sin(a) * sin_comp + np.cos(a) * cos_comp, origin) for a in angles]
+        directions = Vector3Array(np.array([np.sin(a) * sin_comp + np.cos(a) * cos_comp for a in angles]))
+        return list(map(Ray, directions, [origin]*angle_reselution))
     
         
\ No newline at end of file
diff --git a/Raytrace/__init__.py b/Raytrace/__init__.py
index e69de29..9954358 100644
--- a/Raytrace/__init__.py
+++ b/Raytrace/__init__.py
@@ -0,0 +1,5 @@
+from .Meshs import (
+    TriangleMesh as TriangleMesh,
+    BVHMesh as BVHMesh,
+    CompositeMesh as CompositeMesh,
+)
\ No newline at end of file
diff --git a/ShipGen/Generate.py b/ShipGen/Generate.py
index 321ce90..995f14f 100644
--- a/ShipGen/Generate.py
+++ b/ShipGen/Generate.py
@@ -3,9 +3,12 @@
 import numpy as np
 from tqdm import tqdm
 import torch
+from stl import Mesh
+
+from CustomTypes import TriangleArray
+
 from .tools import Guided_Cond_DDPM_Tools as GC_DDPM
 from .tools.HullParameterization import Hull_Parameterization as HP
-from stl import Mesh
 
 data_path = os.path.normpath(os.path.join(__file__,'../data'))
 trained_models_path = os.path.normpath(os.path.join(__file__,'../TrainedModels'))
@@ -77,9 +80,9 @@ def generate_hulls(find_best:int,
             print('Error at hull {}'.format(valid_idx[idxs[i]]))
     return meshs
 
-def mesh_to_triangles(mesh:Mesh) -> np.ndarray:
-    return np.array([mesh.v0, mesh.v1, mesh.v2]).swapaxes(0, 1)
-def mesh_list_to_triangles(mesh_list:list[Mesh]) -> list[np.ndarray]:
+def mesh_to_triangles(mesh:Mesh) -> TriangleArray:
+    return mesh.vectors.view(TriangleArray)
+def mesh_list_to_triangles(mesh_list:list[Mesh]) -> list[TriangleArray]:
     return list(map(mesh_to_triangles, mesh_list))
 
 def normalize_x(mesh:Mesh, scale:float = 10):