diff --git a/examples/01_simple_route.png b/examples/01_simple_route.png
deleted file mode 100644
index 61aa768..0000000
Binary files a/examples/01_simple_route.png and /dev/null differ
diff --git a/examples/02_congestion_resolution.png b/examples/02_congestion_resolution.png
deleted file mode 100644
index 018ae45..0000000
Binary files a/examples/02_congestion_resolution.png and /dev/null differ
diff --git a/examples/03_locked_paths.png b/examples/03_locked_paths.png
index 8995b55..9eb4aec 100644
Binary files a/examples/03_locked_paths.png and b/examples/03_locked_paths.png differ
diff --git a/examples/04_sbends_and_radii.png b/examples/04_sbends_and_radii.png
deleted file mode 100644
index fe3a243..0000000
Binary files a/examples/04_sbends_and_radii.png and /dev/null differ
diff --git a/examples/05_orientation_stress.png b/examples/05_orientation_stress.png
index 035c1ce..7d87622 100644
Binary files a/examples/05_orientation_stress.png and b/examples/05_orientation_stress.png differ
diff --git a/examples/06_bend_collision_models.png b/examples/06_bend_collision_models.png
index 142e346..0dc3212 100644
Binary files a/examples/06_bend_collision_models.png and b/examples/06_bend_collision_models.png differ
diff --git a/examples/07_large_scale_routing.py b/examples/07_large_scale_routing.py
index b20b0da..dc4ad42 100644
--- a/examples/07_large_scale_routing.py
+++ b/examples/07_large_scale_routing.py
@@ -104,7 +104,7 @@ def main() -> None:
 
         # Save plots only for certain iterations to save time
 #        if idx % 20 == 0 or idx == pf.max_iterations - 1:
-        if False:
+        if True:
             # Save a plot of this iteration's result
             fig, ax = plot_routing_results(current_results, obstacles, bounds, netlist=netlist)
             plot_danger_map(danger_map, ax=ax)
diff --git a/examples/08_custom_bend_geometry.png b/examples/08_custom_bend_geometry.png
index 5f51faa..fea9c69 100644
Binary files a/examples/08_custom_bend_geometry.png and b/examples/08_custom_bend_geometry.png differ
diff --git a/examples/09_unroutable_best_effort.png b/examples/09_unroutable_best_effort.png
deleted file mode 100644
index 38c8a2b..0000000
Binary files a/examples/09_unroutable_best_effort.png and /dev/null differ
diff --git a/inire/geometry/collision.py b/inire/geometry/collision.py
index 793d4ec..69b92a9 100644
--- a/inire/geometry/collision.py
+++ b/inire/geometry/collision.py
@@ -96,7 +96,7 @@ class CollisionEngine:
                 f"  Congestion: {m['congestion_tree_queries']} checks\n"
                 f"  Safety Zone: {m['safety_zone_checks']} full intersections performed")
 
-    def add_static_obstacle(self, polygon: Polygon) -> None:
+    def add_static_obstacle(self, polygon: Polygon) -> int:
         obj_id = self._static_id_counter
         self._static_id_counter += 1
         
@@ -115,6 +115,26 @@ class CollisionEngine:
         b = dilated.bounds
         area = (b[2] - b[0]) * (b[3] - b[1])
         self.static_is_rect[obj_id] = (abs(dilated.area - area) < 1e-4)
+        return obj_id
+
+    def remove_static_obstacle(self, obj_id: int) -> None:
+        """
+        Remove a static obstacle by ID.
+        """
+        if obj_id not in self.static_geometries:
+            return
+            
+        bounds = self.static_dilated[obj_id].bounds
+        self.static_index.delete(obj_id, bounds)
+        
+        del self.static_geometries[obj_id]
+        del self.static_dilated[obj_id]
+        del self.static_prepared[obj_id]
+        del self.static_is_rect[obj_id]
+        
+        self.static_tree = None
+        self._static_raw_tree = None
+        self.static_grid = {}
 
     def _ensure_static_tree(self) -> None:
         if self.static_tree is None and self.static_dilated:
@@ -229,26 +249,43 @@ class CollisionEngine:
         tb = result.total_dilated_bounds
         if tb is None: return 0
         self._ensure_dynamic_grid()
-        if not self.dynamic_grid: return 0
+        dynamic_grid = self.dynamic_grid
+        if not dynamic_grid: return 0
         
         cs_inv = self._inv_grid_cell_size
-        gx_min, gy_min = int(tb[0] * cs_inv), int(tb[1] * cs_inv)
-        gx_max, gy_max = int(tb[2] * cs_inv), int(tb[3] * cs_inv)
+        gx_min = int(tb[0] * cs_inv)
+        gy_min = int(tb[1] * cs_inv)
+        gx_max = int(tb[2] * cs_inv)
+        gy_max = int(tb[3] * cs_inv)
         
-        any_possible = False
-        dynamic_grid = self.dynamic_grid
         dynamic_geometries = self.dynamic_geometries
+        
+        # Fast path for single cell
+        if gx_min == gx_max and gy_min == gy_max:
+            cell = (gx_min, gy_min)
+            if cell in dynamic_grid:
+                for obj_id in dynamic_grid[cell]:
+                    if dynamic_geometries[obj_id][0] != net_id:
+                        return self._check_real_congestion(result, net_id)
+            return 0
+
+        # General case
+        any_possible = False
         for gx in range(gx_min, gx_max + 1):
             for gy in range(gy_min, gy_max + 1):
                 cell = (gx, gy)
                 if cell in dynamic_grid:
                     for obj_id in dynamic_grid[cell]:
                         if dynamic_geometries[obj_id][0] != net_id:
-                            any_possible = True; break
+                            any_possible = True
+                            break
                 if any_possible: break
             if any_possible: break
+            
         if not any_possible: return 0
-        
+        return self._check_real_congestion(result, net_id)
+
+    def _check_real_congestion(self, result: ComponentResult, net_id: str) -> int:
         self.metrics['congestion_tree_queries'] += 1
         self._ensure_dynamic_tree()
         if self.dynamic_tree is None: return 0
@@ -259,28 +296,19 @@ class CollisionEngine:
         possible_total = (tb[0] < d_bounds[:, 2]) & (tb[2] > d_bounds[:, 0]) & \
                          (tb[1] < d_bounds[:, 3]) & (tb[3] > d_bounds[:, 1])
         
-        # Filter by net_id (important for negotiated congestion)
         valid_hits = (self._dynamic_net_ids_array != net_id)
         if not numpy.any(possible_total & valid_hits):
              return 0
 
         # 2. Per-polygon AABB check using query on geometries (LAZY triggering)
-        # We only trigger evaluation if total bounds intersect with other nets.
         geoms_to_test = result.dilated_geometry if result.dilated_geometry else result.geometry
         res_indices, tree_indices = self.dynamic_tree.query(geoms_to_test, predicate='intersects')
         
         if tree_indices.size == 0:
              return 0
              
-        # Filter out self-overlaps (from same net)
         hit_net_ids = numpy.take(self._dynamic_net_ids_array, tree_indices)
         valid_geoms_hits = (hit_net_ids != net_id)
-        if not numpy.any(valid_geoms_hits):
-             return 0
-
-        # 3. Real geometry check (Only if AABBs intersect with other nets)
-        # We already have hits from STRtree which are accurate for polygons too.
-        # But wait, query(..., predicate='intersects') ALREADY does real check!
         return int(numpy.sum(valid_geoms_hits))
 
     def _is_in_safety_zone(self, geometry: Polygon, obj_id: int, start_port: Port | None, end_port: Port | None) -> bool:
diff --git a/inire/geometry/components.py b/inire/geometry/components.py
index 742d82c..0808694 100644
--- a/inire/geometry/components.py
+++ b/inire/geometry/components.py
@@ -30,7 +30,7 @@ class ComponentResult:
     __slots__ = (
         '_geometry', '_dilated_geometry', '_proxy_geometry', '_actual_geometry', '_dilated_actual_geometry',
         'end_port', 'length', 'move_type', '_bounds', '_dilated_bounds', 
-        '_total_bounds', '_total_dilated_bounds', '_t_cache', '_total_geom_list', '_offsets', '_coords_cache',
+        '_total_bounds', '_total_dilated_bounds', '_bounds_cached', '_total_geom_list', '_offsets', '_coords_cache',
         '_base_result', '_offset', '_lazy_evaluated', 'rel_gx', 'rel_gy', 'rel_go'
     )
 
@@ -58,11 +58,11 @@ class ComponentResult:
         self.end_port = end_port
         self.length = length
         self.move_type = move_type
-        self._t_cache = {}
         
         self._base_result = _base_result
         self._offset = _offset
         self._lazy_evaluated = False
+        self._bounds_cached = False
         
         if rel_gx is not None:
              self.rel_gx = rel_gx
@@ -138,6 +138,7 @@ class ComponentResult:
                 self._total_bounds = None
                 self._dilated_bounds = None
                 self._total_dilated_bounds = None
+            self._bounds_cached = True
 
     def _ensure_evaluated(self) -> None:
         if self._base_result is None or self._lazy_evaluated:
@@ -145,16 +146,11 @@ class ComponentResult:
         
         # Perform Translation
         dx, dy = self._offset
-        # Base uses its own coords cache
-        base_coords = self._base_result._coords_cache
-        if base_coords is None:
-             self._lazy_evaluated = True
-             return
-             
-        new_coords = base_coords + [dx, dy]
         
-        # Translate ALL geometries at once
-        new_total_arr = shapely.set_coordinates(list(self._base_result._total_geom_list), new_coords)
+        # Use shapely.transform which is vectorized and doesn't modify in-place.
+        # This is MUCH faster than cloning with copy.copy and then set_coordinates.
+        import shapely
+        new_total_arr = shapely.transform(self._base_result._total_geom_list, lambda x: x + [dx, dy])
         new_total = new_total_arr.tolist()
         
         o = self._base_result._offsets
@@ -193,47 +189,53 @@ class ComponentResult:
 
     @property
     def bounds(self) -> numpy.ndarray:
-        if self._bounds is None:
-            if self._base_result is not None:
-                 dx, dy = self._offset
-                 self._bounds = self._base_result.bounds + [dx, dy, dx, dy]
+        if not self._bounds_cached:
+            self._ensure_bounds_evaluated()
         return self._bounds
 
     @property
     def total_bounds(self) -> numpy.ndarray:
-        if self._total_bounds is None:
-            if self._base_result is not None:
-                 dx, dy = self._offset
-                 self._total_bounds = self._base_result.total_bounds + [dx, dy, dx, dy]
+        if not self._bounds_cached:
+            self._ensure_bounds_evaluated()
         return self._total_bounds
 
     @property
     def dilated_bounds(self) -> numpy.ndarray | None:
-        if self._dilated_bounds is None:
-            if self._base_result is not None and self._base_result.dilated_bounds is not None:
-                 dx, dy = self._offset
-                 self._dilated_bounds = self._base_result.dilated_bounds + [dx, dy, dx, dy]
+        if not self._bounds_cached:
+            self._ensure_bounds_evaluated()
         return self._dilated_bounds
 
     @property
     def total_dilated_bounds(self) -> numpy.ndarray | None:
-        if self._total_dilated_bounds is None:
-            if self._base_result is not None and self._base_result.total_dilated_bounds is not None:
-                 dx, dy = self._offset
-                 self._total_dilated_bounds = self._base_result.total_dilated_bounds + [dx, dy, dx, dy]
+        if not self._bounds_cached:
+            self._ensure_bounds_evaluated()
         return self._total_dilated_bounds
 
+    def _ensure_bounds_evaluated(self) -> None:
+        if self._bounds_cached: return
+        base = self._base_result
+        if base is not None:
+             dx, dy = self._offset
+             # Vectorized addition is faster if we avoid creating new lists/arrays in the loop
+             if base._bounds is not None:
+                  self._bounds = base._bounds + [dx, dy, dx, dy]
+             if base._total_bounds is not None:
+                  b = base._total_bounds
+                  self._total_bounds = numpy.array([b[0]+dx, b[1]+dy, b[2]+dx, b[3]+dy])
+             if base._dilated_bounds is not None:
+                  self._dilated_bounds = base._dilated_bounds + [dx, dy, dx, dy]
+             if base._total_dilated_bounds is not None:
+                  b = base._total_dilated_bounds
+                  self._total_dilated_bounds = numpy.array([b[0]+dx, b[1]+dy, b[2]+dx, b[3]+dy])
+        self._bounds_cached = True
+
     def translate(self, dx: float, dy: float, rel_gx: int | None = None, rel_gy: int | None = None, rel_go: int | None = None) -> ComponentResult:
         """
         Create a new ComponentResult translated by (dx, dy).
         """
-        dxr, dyr = round(dx, 3), round(dy, 3)
-        if (dxr, dyr) == (0.0, 0.0):
-            return self
-        if (dxr, dyr) in self._t_cache:
-            return self._t_cache[(dxr, dyr)]
-
-        new_port = Port(self.end_port.x + dx, self.end_port.y + dy, self.end_port.orientation)
+        # Optimized: no internal cache (already cached in router) and no rounding
+        # Also skip snapping since parent and relative move are already snapped
+        new_port = Port(self.end_port.x + dx, self.end_port.y + dy, self.end_port.orientation, snap=False)
         
         # LAZY TRANSLATE
         if self._base_result:
@@ -244,7 +246,7 @@ class ComponentResult:
              base = self
              new_offset = [dx, dy]
 
-        res = ComponentResult(
+        return ComponentResult(
             end_port=new_port,
             length=self.length,
             move_type=self.move_type,
@@ -254,9 +256,6 @@ class ComponentResult:
             rel_gy=rel_gy,
             rel_go=rel_go
         )
-            
-        self._t_cache[(dxr, dyr)] = res
-        return res
 
 
 class Straight:
diff --git a/inire/geometry/primitives.py b/inire/geometry/primitives.py
index ab4ff39..a99e557 100644
--- a/inire/geometry/primitives.py
+++ b/inire/geometry/primitives.py
@@ -25,10 +25,20 @@ class Port:
             x: float,
             y: float,
             orientation: float,
+            snap: bool = True
             ) -> None:
-        self.x = snap_nm(x)
-        self.y = snap_nm(y)
-        self.orientation = float(orientation % 360)
+        if snap:
+            self.x = round(x * 1000) / 1000
+            self.y = round(y * 1000) / 1000
+            # Faster orientation normalization for common cases
+            if 0 <= orientation < 360:
+                self.orientation = float(orientation)
+            else:
+                self.orientation = float(orientation % 360)
+        else:
+            self.x = x
+            self.y = y
+            self.orientation = float(orientation)
 
     def __repr__(self) -> str:
         return f'Port(x={self.x}, y={self.y}, orientation={self.orientation})'
diff --git a/inire/router/astar.py b/inire/router/astar.py
index c0a9019..abc0e84 100644
--- a/inire/router/astar.py
+++ b/inire/router/astar.py
@@ -121,6 +121,7 @@ class AStarRouter:
             return_partial: bool = False,
             store_expanded: bool = False,
             skip_congestion: bool = False,
+            max_cost: float | None = None
             ) -> list[ComponentResult] | None:
         """
         Route a single net using A*.
@@ -155,6 +156,11 @@ class AStarRouter:
 
             current = heapq.heappop(open_set)
             
+            # Cost Pruning (Fail Fast)
+            if max_cost is not None and current.f_cost > max_cost:
+                 self.metrics['pruned_cost'] += 1
+                 continue
+            
             if current.h_cost < best_node.h_cost:
                 best_node = current
 
@@ -177,7 +183,7 @@ class AStarRouter:
                 return self._reconstruct_path(current)
 
             # Expansion
-            self._expand_moves(current, target, net_width, net_id, open_set, closed_set, snap, nodes_expanded, skip_congestion=skip_congestion, inv_snap=inv_snap, parent_state=state)
+            self._expand_moves(current, target, net_width, net_id, open_set, closed_set, snap, nodes_expanded, skip_congestion=skip_congestion, inv_snap=inv_snap, parent_state=state, max_cost=max_cost)
 
         return self._reconstruct_path(best_node) if return_partial else None
 
@@ -193,7 +199,8 @@ class AStarRouter:
             nodes_expanded: int = 0,
             skip_congestion: bool = False,
             inv_snap: float | None = None,
-            parent_state: tuple[int, int, int] | None = None
+            parent_state: tuple[int, int, int] | None = None,
+            max_cost: float | None = None
             ) -> None:
         cp = current.port
         if inv_snap is None: inv_snap = 1.0 / snap
@@ -214,7 +221,7 @@ class AStarRouter:
         if proj_t > 0 and abs(perp_t) < 1e-3 and abs(cp.orientation - target.orientation) < 0.1:
              max_reach = self.cost_evaluator.collision_engine.ray_cast(cp, cp.orientation, proj_t + 1.0)
              if max_reach >= proj_t - 0.01:
-                  self._process_move(current, target, net_width, net_id, open_set, closed_set, snap, f'S{proj_t}', 'S', (proj_t,), skip_congestion, inv_snap=inv_snap, snap_to_grid=False, parent_state=parent_state)
+                  self._process_move(current, target, net_width, net_id, open_set, closed_set, snap, f'S{proj_t}', 'S', (proj_t,), skip_congestion, inv_snap=inv_snap, snap_to_grid=False, parent_state=parent_state, max_cost=max_cost)
 
         # 2. VISIBILITY JUMPS & MAX REACH
         max_reach = self.cost_evaluator.collision_engine.ray_cast(cp, cp.orientation, self.config.max_straight_length)
@@ -261,7 +268,7 @@ class AStarRouter:
                             if s_l <= max_reach and s_l > 0.1: straight_lengths.add(s_l)
 
         for length in sorted(straight_lengths, reverse=True):
-            self._process_move(current, target, net_width, net_id, open_set, closed_set, snap, f'S{length}', 'S', (length,), skip_congestion, inv_snap=inv_snap, parent_state=parent_state)
+            self._process_move(current, target, net_width, net_id, open_set, closed_set, snap, f'S{length}', 'S', (length,), skip_congestion, inv_snap=inv_snap, parent_state=parent_state, max_cost=max_cost)
 
         # 3. BENDS & SBENDS
         angle_to_target = numpy.degrees(numpy.arctan2(target.y - cp.y, target.x - cp.x))
@@ -275,7 +282,7 @@ class AStarRouter:
                     new_diff = (angle_to_target - new_ori + 180) % 360 - 180
                     if abs(new_diff) > 135:
                         continue
-                self._process_move(current, target, net_width, net_id, open_set, closed_set, snap, f'B{radius}{direction}', 'B', (radius, direction), skip_congestion, inv_snap=inv_snap, parent_state=parent_state)
+                self._process_move(current, target, net_width, net_id, open_set, closed_set, snap, f'B{radius}{direction}', 'B', (radius, direction), skip_congestion, inv_snap=inv_snap, parent_state=parent_state, max_cost=max_cost)
 
         # 4. SBENDS
         max_sbend_r = max(self.config.sbend_radii) if self.config.sbend_radii else 0
@@ -298,7 +305,7 @@ class AStarRouter:
             for offset in sorted(offsets):
                 for radius in self.config.sbend_radii:
                     if abs(offset) >= 2 * radius: continue
-                    self._process_move(current, target, net_width, net_id, open_set, closed_set, snap, f'SB{offset}R{radius}', 'SB', (offset, radius), skip_congestion, inv_snap=inv_snap, parent_state=parent_state)
+                    self._process_move(current, target, net_width, net_id, open_set, closed_set, snap, f'SB{offset}R{radius}', 'SB', (offset, radius), skip_congestion, inv_snap=inv_snap, parent_state=parent_state, max_cost=max_cost)
 
     def _process_move(
             self,
@@ -315,7 +322,8 @@ class AStarRouter:
             skip_congestion: bool,
             inv_snap: float | None = None,
             snap_to_grid: bool = True,
-            parent_state: tuple[int, int, int] | None = None
+            parent_state: tuple[int, int, int] | None = None,
+            max_cost: float | None = None
             ) -> None:
         cp = parent.port
         if inv_snap is None: inv_snap = 1.0 / snap
@@ -333,7 +341,7 @@ class AStarRouter:
         if abs_key in self._move_cache:
             res = self._move_cache[abs_key]
             move_radius = params[0] if move_class == 'B' else (params[1] if move_class == 'SB' else None)
-            self._add_node(parent, res, target, net_width, net_id, open_set, closed_set, move_type, move_radius=move_radius, snap=snap, skip_congestion=skip_congestion, inv_snap=inv_snap, parent_state=parent_state)
+            self._add_node(parent, res, target, net_width, net_id, open_set, closed_set, move_type, move_radius=move_radius, snap=snap, skip_congestion=skip_congestion, inv_snap=inv_snap, parent_state=parent_state, max_cost=max_cost)
             return
 
         rel_key = (base_ori, move_class, params, net_width, self.config.bend_collision_type, self._self_dilation, snap_to_grid)
@@ -344,7 +352,6 @@ class AStarRouter:
             
         if rel_key in self._move_cache:
             res_rel = self._move_cache[rel_key]
-            res = res_rel.translate(cp.x, cp.y, rel_gx=res_rel.rel_gx + gx, rel_gy=res_rel.rel_gy + gy, rel_go=res_rel.rel_go)
         else:
             try:
                 p0 = Port(0, 0, base_ori)
@@ -357,13 +364,13 @@ class AStarRouter:
                 else:
                     return
                 self._move_cache[rel_key] = res_rel
-                res = res_rel.translate(cp.x, cp.y, rel_gx=res_rel.rel_gx + gx, rel_gy=res_rel.rel_gy + gy, rel_go=res_rel.rel_go)
             except (ValueError, ZeroDivisionError):
                 return
 
+        res = res_rel.translate(cp.x, cp.y, rel_gx=res_rel.rel_gx + gx, rel_gy=res_rel.rel_gy + gy, rel_go=res_rel.rel_go)
         self._move_cache[abs_key] = res
         move_radius = params[0] if move_class == 'B' else (params[1] if move_class == 'SB' else None)
-        self._add_node(parent, res, target, net_width, net_id, open_set, closed_set, move_type, move_radius=move_radius, snap=snap, skip_congestion=skip_congestion, inv_snap=inv_snap, parent_state=parent_state)
+        self._add_node(parent, res, target, net_width, net_id, open_set, closed_set, move_type, move_radius=move_radius, snap=snap, skip_congestion=skip_congestion, inv_snap=inv_snap, parent_state=parent_state, max_cost=max_cost)
 
     def _add_node(
             self,
@@ -379,7 +386,8 @@ class AStarRouter:
             snap: float = 1.0,
             skip_congestion: bool = False,
             inv_snap: float | None = None,
-            parent_state: tuple[int, int, int] | None = None
+            parent_state: tuple[int, int, int] | None = None,
+            max_cost: float | None = None
             ) -> None:
         self.metrics['moves_generated'] += 1
         state = (result.rel_gx, result.rel_gy, result.rel_go)
@@ -400,6 +408,15 @@ class AStarRouter:
             self.metrics['pruned_hard_collision'] += 1
             return
 
+        new_g_cost = parent.g_cost + result.length
+        
+        # Pre-check cost pruning before evaluation (using heuristic)
+        if max_cost is not None:
+             new_h_cost = self.cost_evaluator.h_manhattan(end_p, target)
+             if new_g_cost + new_h_cost > max_cost:
+                  self.metrics['pruned_cost'] += 1
+                  return
+
         is_static_safe = (cache_key in self._static_safe_cache)
         if not is_static_safe:
             ce = self.cost_evaluator.collision_engine
diff --git a/inire/router/cost.py b/inire/router/cost.py
index 165cce9..e6d45c4 100644
--- a/inire/router/cost.py
+++ b/inire/router/cost.py
@@ -106,16 +106,14 @@ class CostEvaluator:
         """
         Heuristic: weighted Manhattan distance + mandatory turn penalties.
         """
-        tx = target.x
-        ty = target.y
+        tx, ty = target.x, target.y
         t_ori = target.orientation
-        t_cos = self._target_cos
-        t_sin = self._target_sin
         
+        # Avoid repeated trig for target orientation
         if abs(tx - self._target_x) > 1e-6 or abs(ty - self._target_y) > 1e-6:
-             rad = np.radians(t_ori)
-             t_cos = np.cos(rad)
-             t_sin = np.sin(rad)
+             self.set_target(target)
+        
+        t_cos, t_sin = self._target_cos, self._target_sin
 
         dx = abs(current.x - tx)
         dy = abs(current.y - ty)
@@ -125,7 +123,9 @@ class CostEvaluator:
         penalty = 0.0
 
         # 1. Orientation Difference
-        diff = abs(current.orientation - t_ori) % 360
+        # Optimization: use integer comparison for common orientations
+        curr_ori = current.orientation
+        diff = abs(curr_ori - t_ori) % 360
         if diff > 0.1:
             if abs(diff - 180) < 0.1:
                 penalty += 2 * bp
@@ -143,8 +143,16 @@ class CostEvaluator:
             penalty += 2 * bp
 
         # 3. Traveling Away
-        curr_rad = np.radians(current.orientation)
-        move_proj = v_dx * np.cos(curr_rad) + v_dy * np.sin(curr_rad)
+        # Optimization: avoid np.radians/cos/sin if current_ori is standard 0,90,180,270
+        if curr_ori == 0: c_cos, c_sin = 1.0, 0.0
+        elif curr_ori == 90: c_cos, c_sin = 0.0, 1.0
+        elif curr_ori == 180: c_cos, c_sin = -1.0, 0.0
+        elif curr_ori == 270: c_cos, c_sin = 0.0, -1.0
+        else:
+            curr_rad = np.radians(curr_ori)
+            c_cos, c_sin = np.cos(curr_rad), np.sin(curr_rad)
+            
+        move_proj = v_dx * c_cos + v_dy * c_sin
         if move_proj < -0.1:
             penalty += 2 * bp
 
diff --git a/inire/router/pathfinder.py b/inire/router/pathfinder.py
index 94fb77d..107ea7c 100644
--- a/inire/router/pathfinder.py
+++ b/inire/router/pathfinder.py
@@ -4,7 +4,7 @@ import logging
 import time
 import random
 from dataclasses import dataclass
-from typing import TYPE_CHECKING, Callable
+from typing import TYPE_CHECKING, Callable, Literal, Any
 
 if TYPE_CHECKING:
     from inire.geometry.components import ComponentResult
@@ -40,7 +40,7 @@ class PathFinder:
     """
     Multi-net router using Negotiated Congestion.
     """
-    __slots__ = ('router', 'cost_evaluator', 'max_iterations', 'base_congestion_penalty', 'use_tiered_strategy', 'congestion_multiplier', 'accumulated_expanded_nodes')
+    __slots__ = ('router', 'cost_evaluator', 'max_iterations', 'base_congestion_penalty', 'use_tiered_strategy', 'congestion_multiplier', 'accumulated_expanded_nodes', 'warm_start')
 
     router: AStarRouter
     """ The A* search engine """
@@ -60,6 +60,9 @@ class PathFinder:
     use_tiered_strategy: bool
     """ If True, use simpler collision models in early iterations for speed """
 
+    warm_start: Literal['shortest', 'longest', 'user'] | None
+    """ Heuristic sorting for the initial greedy pass """
+
     def __init__(
             self,
             router: AStarRouter,
@@ -68,6 +71,7 @@ class PathFinder:
             base_congestion_penalty: float = 100.0,
             congestion_multiplier: float = 1.5,
             use_tiered_strategy: bool = True,
+            warm_start: Literal['shortest', 'longest', 'user'] | None = 'shortest',
             ) -> None:
         """
         Initialize the PathFinder.
@@ -79,6 +83,7 @@ class PathFinder:
             base_congestion_penalty: Starting penalty for overlaps.
             congestion_multiplier: Multiplier for congestion penalty per iteration.
             use_tiered_strategy: Whether to use simplified collision models in early iterations.
+            warm_start: Initial ordering strategy for a fast greedy pass.
         """
         self.router = router
         self.cost_evaluator = cost_evaluator
@@ -86,8 +91,59 @@ class PathFinder:
         self.base_congestion_penalty = base_congestion_penalty
         self.congestion_multiplier = congestion_multiplier
         self.use_tiered_strategy = use_tiered_strategy
+        self.warm_start = warm_start
         self.accumulated_expanded_nodes: list[tuple[float, float, float]] = []
 
+    def _perform_greedy_pass(
+            self,
+            netlist: dict[str, tuple[Port, Port]],
+            net_widths: dict[str, float],
+            order: Literal['shortest', 'longest', 'user']
+            ) -> dict[str, list[ComponentResult]]:
+        """
+        Internal greedy pass: route nets sequentially and freeze them as static.
+        """
+        all_net_ids = list(netlist.keys())
+        if order != 'user':
+            def get_dist(nid):
+                s, t = netlist[nid]
+                return abs(t.x - s.x) + abs(t.y - s.y)
+            all_net_ids.sort(key=get_dist, reverse=(order == 'longest'))
+            
+        greedy_paths = {}
+        temp_obj_ids = []
+        
+        logger.info(f"PathFinder: Starting Greedy Warm-Start ({order} order)...")
+        
+        for net_id in all_net_ids:
+            start, target = netlist[net_id]
+            width = net_widths.get(net_id, 2.0)
+            
+            # Heuristic max cost for fail-fast
+            h_start = self.cost_evaluator.h_manhattan(start, target)
+            max_cost_limit = max(h_start * 3.0, 2000.0) 
+            
+            path = self.router.route(
+                start, target, width, net_id=net_id, 
+                skip_congestion=True, max_cost=max_cost_limit
+            )
+            
+            if path:
+                greedy_paths[net_id] = path
+                # Freeze as static
+                for res in path:
+                    geoms = res.actual_geometry if res.actual_geometry is not None else res.geometry
+                    for poly in geoms:
+                        obj_id = self.cost_evaluator.collision_engine.add_static_obstacle(poly)
+                        temp_obj_ids.append(obj_id)
+        
+        # Clean up temporary static obstacles
+        for obj_id in temp_obj_ids:
+            self.cost_evaluator.collision_engine.remove_static_obstacle(obj_id)
+            
+        logger.info(f"PathFinder: Greedy Warm-Start finished. Seeding {len(greedy_paths)}/{len(netlist)} nets.")
+        return greedy_paths
+
     def route_all(
             self,
             netlist: dict[str, tuple[Port, Port]],
@@ -95,6 +151,8 @@ class PathFinder:
             store_expanded: bool = False,
             iteration_callback: Callable[[int, dict[str, RoutingResult]], None] | None = None,
             shuffle_nets: bool = False,
+            sort_nets: Literal['shortest', 'longest', 'user', None] = None,
+            initial_paths: dict[str, list[ComponentResult]] | None = None,
             seed: int | None = None,
             ) -> dict[str, RoutingResult]:
         """
@@ -106,6 +164,8 @@ class PathFinder:
             store_expanded: Whether to store expanded nodes for ALL iterations and nets.
             iteration_callback: Optional callback(iteration_idx, current_results).
             shuffle_nets: Whether to randomize the order of nets each iteration.
+            sort_nets: Heuristic sorting for the initial iteration order (overrides self.warm_start).
+            initial_paths: Pre-computed paths to use for Iteration 0 (overrides warm_start).
             seed: Optional seed for randomization (enables reproducibility).
 
         Returns:
@@ -120,6 +180,21 @@ class PathFinder:
         session_timeout = max(60.0, 10.0 * num_nets * self.max_iterations)
 
         all_net_ids = list(netlist.keys())
+        
+        # Determine initial paths (Warm Start)
+        if initial_paths is None:
+             ws_order = sort_nets if sort_nets is not None else self.warm_start
+             if ws_order is not None:
+                  initial_paths = self._perform_greedy_pass(netlist, net_widths, ws_order)
+
+        # Apply initial sorting heuristic if requested (for the main NC loop)
+        if sort_nets:
+            def get_dist(nid):
+                s, t = netlist[nid]
+                return abs(t.x - s.x) + abs(t.y - s.y)
+            
+            if sort_nets != 'user':
+                 all_net_ids.sort(key=get_dist, reverse=(sort_nets == 'longest'))
 
         for iteration in range(self.max_iterations):
             any_congestion = False
@@ -148,36 +223,39 @@ class PathFinder:
                 # 1. Rip-up existing path
                 self.cost_evaluator.collision_engine.remove_path(net_id)
 
-                # 2. Reroute with current congestion info
-                target_coll_model = self.router.config.bend_collision_type
-                coll_model = target_coll_model
-                skip_cong = False
-                if self.use_tiered_strategy and iteration == 0:
-                    skip_cong = True
-                    if target_coll_model == "arc":
-                        coll_model = "clipped_bbox"
+                # 2. Reroute or Use Initial Path
+                path = None
                 
-                # Dynamic node limit: increase if it failed previously
-                base_node_limit = self.router.config.node_limit
-                current_node_limit = base_node_limit
-                if net_id in results and not results[net_id].reached_target:
-                     current_node_limit = base_node_limit * (iteration + 1)
-                
-                net_start = time.monotonic()
-                
-                # Temporarily override node_limit
-                original_limit = self.router.node_limit
-                self.router.node_limit = current_node_limit
-                
-                path = self.router.route(start, target, width, net_id=net_id, bend_collision_type=coll_model, return_partial=True, store_expanded=store_expanded, skip_congestion=skip_cong)
-                
-                if store_expanded and self.router.last_expanded_nodes:
-                    self.accumulated_expanded_nodes.extend(self.router.last_expanded_nodes)
+                # Warm Start Logic: Use provided path for Iteration 0
+                if iteration == 0 and initial_paths and net_id in initial_paths:
+                    path = initial_paths[net_id]
+                    logger.debug(f'  Net {net_id} used Warm Start path.')
+                else:
+                    # Standard Routing Logic
+                    target_coll_model = self.router.config.bend_collision_type
+                    coll_model = target_coll_model
+                    skip_cong = False
+                    if self.use_tiered_strategy and iteration == 0:
+                        skip_cong = True
+                        if target_coll_model == "arc":
+                            coll_model = "clipped_bbox"
+                    
+                    base_node_limit = self.router.config.node_limit
+                    current_node_limit = base_node_limit
+                    if net_id in results and not results[net_id].reached_target:
+                         current_node_limit = base_node_limit * (iteration + 1)
+                    
+                    net_start = time.monotonic()
+                    original_limit = self.router.node_limit
+                    self.router.node_limit = current_node_limit
+                    
+                    path = self.router.route(start, target, width, net_id=net_id, bend_collision_type=coll_model, return_partial=True, store_expanded=store_expanded, skip_congestion=skip_cong)
+                    
+                    if store_expanded and self.router.last_expanded_nodes:
+                        self.accumulated_expanded_nodes.extend(self.router.last_expanded_nodes)
 
-                # Restore
-                self.router.node_limit = original_limit
-                
-                logger.debug(f'  Net {net_id} routed in {time.monotonic() - net_start:.4f}s using {coll_model}')
+                    self.router.node_limit = original_limit
+                    logger.debug(f'  Net {net_id} routed in {time.monotonic() - net_start:.4f}s using {coll_model}')
 
                 if path:
                     # Check if reached exactly