diff --git a/benchmarks/benchmarks/common.py b/benchmarks/benchmarks/common.py
index 06207e47d..e2e1dfde8 100644
--- a/benchmarks/benchmarks/common.py
+++ b/benchmarks/benchmarks/common.py
@@ -58,6 +58,13 @@ def get_xr_dataarray(
     elif type == "dask":
         import dask.array as da
         z = da.from_array(z, chunks=(max(1, ny // 2), max(1, nx // 2)))
+    elif type == "dask+cupy":
+        if not has_cuda_and_cupy():
+            raise NotImplementedError()
+        import cupy
+        import dask.array as da
+        z = da.from_array(
+            cupy.asarray(z), chunks=(max(1, ny // 2), max(1, nx // 2)))
     else:
         raise RuntimeError(f"Unrecognised type {type}")
 
diff --git a/benchmarks/benchmarks/pathfinding.py b/benchmarks/benchmarks/pathfinding.py
index 78ca38dd3..e1ecaf551 100644
--- a/benchmarks/benchmarks/pathfinding.py
+++ b/benchmarks/benchmarks/pathfinding.py
@@ -2,6 +2,7 @@
 import xarray as xr
 
 from xrspatial.pathfinding import a_star_search, multi_stop_search
+from xrspatial.utils import has_cuda_and_cupy
 
 from .common import get_xr_dataarray
 
@@ -36,12 +37,12 @@ def peakmem_multi_stop_search(self):
 
 
 class AStarSearch:
-    params = ([100, 300, 1000], [4, 8], ["numpy", "cupy", "dask"])
+    params = ([100, 300, 1000], [4, 8], ["numpy", "cupy", "dask", "dask+cupy"])
     param_names = ("nx", "connectivity", "type")
 
     def setup(self, nx, connectivity, type):
-        if type == "dask" and nx > 300:
-            # The dask backend is a pure-Python sparse A* that loads
+        if type in ("dask", "dask+cupy") and nx > 300:
+            # The dask backends run a pure-Python sparse A* that loads
             # chunks on demand; at nx=1000 a single call takes ~4 s,
             # which would dominate the suite's runtime.
             raise NotImplementedError()
@@ -50,7 +51,7 @@ def setup(self, nx, connectivity, type):
         self.start = self.agg.y[0], self.agg.x[0]
         self.goal = self.agg.y[-1], self.agg.x[-1]
         # snap_start/snap_goal raise on dask-backed arrays by design
-        self.snap = type != "dask"
+        self.snap = type not in ("dask", "dask+cupy")
 
     def time_a_star_search(self, nx, connectivity, type):
         a_star_search(
@@ -60,6 +61,64 @@ def time_a_star_search(self, nx, connectivity, type):
         )
 
 
+class AStarSearchObstacles:
+    """A* through a barrier field with friction (issue #3705).
+
+    The open-grid AStarSearch benchmark is A*'s best case: the heuristic
+    walks almost straight to the goal.  Walls with alternating gaps force
+    a serpentine route, so the frontier and visited set actually grow,
+    and the friction surface exercises the friction-weighted cost path.
+    Regressions in barrier handling or frontier bookkeeping are invisible
+    to the open-grid benchmark but show up here.
+    """
+
+    params = ([100, 300], ["numpy", "cupy", "dask"])
+    param_names = ("nx", "type")
+
+    def setup(self, nx, type):
+        ny = nx // 2
+        z = np.ones((ny, nx), dtype=np.float64)
+        # Vertical walls (value 0) with a one-cell gap alternating
+        # between the top and bottom row
+        step = max(2, nx // 10)
+        for wi, col in enumerate(range(step, nx - 1, step)):
+            z[:, col] = 0.0
+            if wi % 2 == 0:
+                z[0, col] = 1.0
+            else:
+                z[-1, col] = 1.0
+        # Left-to-right friction gradient
+        f = 1.0 + np.linspace(0.0, 4.0, nx)[np.newaxis, :] * np.ones(
+            (ny, 1), dtype=np.float64)
+
+        if type == "cupy":
+            if not has_cuda_and_cupy():
+                raise NotImplementedError()
+            import cupy
+            z = cupy.asarray(z)
+            f = cupy.asarray(f)
+        elif type == "dask":
+            import dask.array as da
+            chunks = (max(1, ny // 2), max(1, nx // 2))
+            z = da.from_array(z, chunks=chunks)
+            f = da.from_array(f, chunks=chunks)
+
+        y = np.linspace(ny - 1, 0, ny)
+        x = np.linspace(0, nx - 1, nx)
+        self.agg = xr.DataArray(z, coords=dict(y=y, x=x), dims=["y", "x"])
+        self.friction = xr.DataArray(
+            f, coords=dict(y=y, x=x), dims=["y", "x"])
+        self.start = float(y[-1]), float(x[0])
+        self.goal = float(y[0]), float(x[-1])
+
+    def time_a_star_search_barriers(self, nx, type):
+        a_star_search(self.agg, self.start, self.goal, barriers=[0])
+
+    def time_a_star_search_barriers_friction(self, nx, type):
+        a_star_search(self.agg, self.start, self.goal, barriers=[0],
+                      friction=self.friction)
+
+
 class MultiStopSearch:
     params = ([100, 300], ["numpy", "cupy", "dask"])
     param_names = ("nx", "type")
diff --git a/xrspatial/pathfinding.py b/xrspatial/pathfinding.py
index f5662750a..4a077f2ed 100644
--- a/xrspatial/pathfinding.py
+++ b/xrspatial/pathfinding.py
@@ -35,6 +35,11 @@ def _validate_surface_dims(surface, x, y, func_name):
             f"got {surface.dims}. Pass the actual dimension names via the "
             f"`x=` and `y=` parameters."
         )
+    if 0 in surface.shape:
+        raise ValueError(
+            f"{func_name}(): `surface` is empty (shape {surface.shape}); "
+            f"pathfinding requires at least one cell."
+        )
 
 
 def _validate_barriers(barriers):
@@ -283,8 +288,8 @@ def _a_star_search(data, path_img, start_py, start_px, goal_py, goal_px,
 
     # parent of the (i, j) pixel is the pixel at
     # (parent_ys[i, j], parent_xs[i, j])
-    parent_ys = np.ones((height, width), dtype=np.int64) * NONE
-    parent_xs = np.ones((height, width), dtype=np.int64) * NONE
+    parent_ys = np.full((height, width), NONE, dtype=np.int64)
+    parent_xs = np.full((height, width), NONE, dtype=np.int64)
 
     # parent of start is itself
     parent_ys[start_py, start_px] = start_py
@@ -651,6 +656,28 @@ def get(self, key, loader):
         return value
 
 
+# Minimum-positive-friction values, keyed by the dask array's graph token
+# (content-derived, so a mutated or different raster gets a fresh entry).
+# multi_stop_search routes N-1 segments through a_star_search with the same
+# friction raster; without this, every segment re-runs a full nanmin scan
+# (#3660 follow-up, bundled in #3705).
+_DASK_FMIN_CACHE_SIZE = 16
+_dask_fmin_cache = OrderedDict()
+
+
+def _dask_friction_fmin(friction_data):
+    key = friction_data.name
+    if key in _dask_fmin_cache:
+        _dask_fmin_cache.move_to_end(key)
+        return _dask_fmin_cache[key]
+    positive_friction = da.where(friction_data > 0, friction_data, np.inf)
+    f_min = float(da.nanmin(positive_friction).compute())
+    if len(_dask_fmin_cache) >= _DASK_FMIN_CACHE_SIZE:
+        _dask_fmin_cache.popitem(last=False)
+    _dask_fmin_cache[key] = f_min
+    return f_min
+
+
 # ---------------------------------------------------------------------------
 # Sparse dask A*
 # ---------------------------------------------------------------------------
@@ -738,6 +765,11 @@ def _get_value(da_arr, cache, r, c):
 
         g_u = g_cost[(py, px)]
 
+        # Friction at the popped node is invariant across its neighbors;
+        # fetch it once instead of once per neighbor.
+        if use_friction:
+            f_u_val = _get_value(friction_da, friction_cache, py, px)
+
         for i in range(n_neighbors):
             ny = py + int(dy[i])
             nx = px + int(dx[i])
@@ -752,7 +784,6 @@ def _get_value(da_arr, cache, r, c):
                 continue
 
             if use_friction:
-                f_u_val = _get_value(friction_da, friction_cache, py, px)
                 f_v_val = _get_value(friction_da, friction_cache, ny, nx)
                 if not (np.isfinite(f_v_val) and f_v_val > 0.0):
                     continue
@@ -1065,12 +1096,15 @@ def a_star_search(surface: xr.DataArray,
                 "snap_goal is not supported with dask-backed arrays; "
                 "ensure the goal pixel is valid before calling a_star_search"
             )
-        # Single-pixel crossability check via .compute()
-        start_val = float(surface_data[start_py, start_px].compute())
-        if _is_not_crossable_py(start_val, barriers):
+        # Single-pixel crossability checks, batched into one compute so
+        # the scheduler round-trip happens once instead of twice.
+        start_val, goal_val = dask.compute(
+            surface_data[start_py, start_px],
+            surface_data[goal_py, goal_px],
+        )
+        if _is_not_crossable_py(float(start_val), barriers):
             warnings.warn("Start at a non crossable location", Warning)
-        goal_val = float(surface_data[goal_py, goal_px].compute())
-        if _is_not_crossable_py(goal_val, barriers):
+        if _is_not_crossable_py(float(goal_val), barriers):
             warnings.warn("End at a non crossable location", Warning)
     elif _is_cupy_backend:
         # CuPy: use .get() for scalar access in numpy-land
@@ -1124,10 +1158,9 @@ def a_star_search(surface: xr.DataArray,
             else:
                 friction_data = da.from_array(
                     friction_data, chunks=surface_data.chunks)
-            # f_min via dask (same pattern as cost_distance)
-            positive_friction = da.where(
-                friction_data > 0, friction_data, np.inf)
-            f_min = float(da.nanmin(positive_friction).compute())
+            # f_min via dask (same pattern as cost_distance), cached so
+            # multi-segment routes scan the friction raster only once
+            f_min = _dask_friction_fmin(friction_data)
             if not (np.isfinite(f_min) and f_min > 0):
                 raise ValueError("friction has no positive finite values")
         else:
@@ -1356,6 +1389,10 @@ def _held_karp(dist, start, end):
 def _nearest_neighbor_2opt(dist, start, end):
     """Heuristic TSP for large N: nearest-neighbor + 2-opt with fixed endpoints.
 
+    A tour containing an unavoidable inf edge (an unreachable waypoint
+    pair) still receives finite local improvements; its total cost stays
+    inf.
+
     Parameters
     ----------
     dist : 2-D array-like, shape (N, N)
@@ -1380,22 +1417,56 @@ def _nearest_neighbor_2opt(dist, start, end):
         cur = nearest
     tour.append(end)
 
-    # 2-opt local search (only swap interior segment, keep endpoints fixed)
-    def _tour_cost(t):
-        return sum(dist[t[i]][t[i + 1]] for i in range(len(t) - 1))
+    # 2-opt local search (only swap interior segment, keep endpoints fixed).
+    # Candidate moves are scored in O(1) from prefix sums over the current
+    # tour's forward and backward edge costs, instead of re-summing the whole
+    # tour per candidate.  Reversing tour[i:j+1] removes the forward edges
+    # i-1..j and adds the two new boundary edges plus the interior edges
+    # traversed backward (which differ from the forward ones when *dist* is
+    # asymmetric, e.g. with snapped waypoints).  Infinite edges (unreachable
+    # waypoint pairs) are tracked as counts next to the sums of finite edges,
+    # because subtracting prefix sums that both contain inf gives nan and
+    # would silently reject valid moves.
+    INF = float('inf')
 
+    def _prefix_sums(t):
+        m = len(t)
+        fwd_sum = [0.0] * m
+        fwd_inf = [0] * m
+        bwd_sum = [0.0] * m
+        bwd_inf = [0] * m
+        for k in range(m - 1):
+            fe = dist[t[k]][t[k + 1]]
+            be = dist[t[k + 1]][t[k]]
+            fwd_sum[k + 1] = fwd_sum[k] + (0.0 if fe == INF else fe)
+            fwd_inf[k + 1] = fwd_inf[k] + (fe == INF)
+            bwd_sum[k + 1] = bwd_sum[k] + (0.0 if be == INF else be)
+            bwd_inf[k + 1] = bwd_inf[k] + (be == INF)
+        return fwd_sum, fwd_inf, bwd_sum, bwd_inf
+
+    def _edge_range(sums, infs, lo, hi):
+        # Sum of tour edges lo..hi-1 (INF if the range has an inf edge)
+        if infs[hi] - infs[lo]:
+            return INF
+        return sums[hi] - sums[lo]
+
+    fwd_sum, fwd_inf, bwd_sum, bwd_inf = _prefix_sums(tour)
     improved = True
     while improved:
         improved = False
         for i in range(1, len(tour) - 2):
             for j in range(i + 1, len(tour) - 1):
-                # Reverse segment tour[i:j+1]
-                new_tour = tour[:i] + tour[i:j + 1][::-1] + tour[j + 1:]
-                if _tour_cost(new_tour) < _tour_cost(tour):
-                    tour = new_tour
+                removed = _edge_range(fwd_sum, fwd_inf, i - 1, j + 1)
+                added = (dist[tour[i - 1]][tour[j]]
+                         + _edge_range(bwd_sum, bwd_inf, i, j)
+                         + dist[tour[i]][tour[j + 1]])
+                if added < removed:
+                    # Reverse segment tour[i:j+1]
+                    tour = tour[:i] + tour[i:j + 1][::-1] + tour[j + 1:]
+                    fwd_sum, fwd_inf, bwd_sum, bwd_inf = _prefix_sums(tour)
                     improved = True
 
-    return tour, _tour_cost(tour)
+    return tour, _edge_range(fwd_sum, fwd_inf, 0, len(tour) - 1)
 
 
 def _segment_to_numpy(seg_data):
diff --git a/xrspatial/tests/test_pathfinding.py b/xrspatial/tests/test_pathfinding.py
index 1ce6d27bb..b645b17ea 100644
--- a/xrspatial/tests/test_pathfinding.py
+++ b/xrspatial/tests/test_pathfinding.py
@@ -2312,3 +2312,205 @@ def test_barriers_default_none_matches_empty_list():
     m_none = multi_stop_search(agg, [start, goal], barriers=None)
     np.testing.assert_allclose(
         m_none.values, r_default.values, equal_nan=True)
+
+
+# ---------------------------------------------------------------------------
+# LOW-backlog cleanup tests (#3705)
+# ---------------------------------------------------------------------------
+
+def test_chunk_cache_eviction():
+    """_ChunkCache evicts the least-recently-used entry at maxsize."""
+    from collections import Counter
+
+    from xrspatial.pathfinding import _ChunkCache
+
+    loads = Counter()
+
+    def make_loader(key):
+        def load():
+            loads[key] += 1
+            return key.upper()
+        return load
+
+    cache = _ChunkCache(maxsize=2)
+    assert cache.get('a', make_loader('a')) == 'A'
+    assert cache.get('b', make_loader('b')) == 'B'
+    # Hit on 'a' refreshes its LRU position without reloading
+    assert cache.get('a', make_loader('a')) == 'A'
+    assert loads['a'] == 1
+
+    # Third key evicts 'b' (least recently used), not 'a'
+    assert cache.get('c', make_loader('c')) == 'C'
+    assert cache.get('a', make_loader('a')) == 'A'
+    assert loads['a'] == 1
+    assert cache.get('b', make_loader('b')) == 'B'
+    assert loads['b'] == 2
+
+
+def test_held_karp_two_cities():
+    """n == 2 shortcut returns the direct edge in both directions."""
+    dist = [[0, 7], [3, 0]]
+    order, cost = _held_karp(dist, 0, 1)
+    assert order == [0, 1]
+    assert cost == 7
+    order, cost = _held_karp(dist, 1, 0)
+    assert order == [1, 0]
+    assert cost == 3
+
+
+def test_nearest_neighbor_2opt_improvement_branch():
+    """A tour the greedy construction gets wrong is fixed by 2-opt.
+
+    Nearest-neighbor from 0 picks 1 (cost 1) then 2 then 3, giving
+    [0, 1, 2, 3] with cost 1 + 1 + 10 = 12.  Reversing the interior
+    gives [0, 2, 1, 3] with cost 2 + 1 + 1 = 4, so the improvement
+    branch must fire.
+    """
+    dist = [
+        [0, 1, 2, 50],
+        [1, 0, 1, 1],
+        [2, 1, 0, 10],
+        [50, 1, 10, 0],
+    ]
+    order, cost = _nearest_neighbor_2opt(dist, 0, 3)
+    assert order == [0, 2, 1, 3]
+    assert cost == 4
+
+
+def test_nearest_neighbor_2opt_asymmetric_consistency():
+    """Delta scoring stays exact on asymmetric matrices.
+
+    The returned cost must equal the cost recomputed from the returned
+    order, and never exceed the greedy nearest-neighbor tour it started
+    from.
+    """
+    def tour_cost(dist, t):
+        return sum(dist[t[k]][t[k + 1]] for k in range(len(t) - 1))
+
+    def greedy(dist, start, end):
+        n = len(dist)
+        remaining = set(range(n)) - {start, end}
+        tour, cur = [start], start
+        while remaining:
+            cur = min(remaining, key=lambda c: dist[cur][c])
+            tour.append(cur)
+            remaining.remove(cur)
+        return tour + [end]
+
+    rng = np.random.default_rng(3705)
+    for n in (5, 8):
+        for _ in range(5):
+            dist = rng.integers(1, 20, size=(n, n)).astype(float).tolist()
+            for i in range(n):
+                dist[i][i] = 0.0
+            order, cost = _nearest_neighbor_2opt(dist, 0, n - 1)
+            assert order[0] == 0 and order[-1] == n - 1
+            assert sorted(order) == list(range(n))
+            assert cost == pytest.approx(tour_cost(dist, order))
+            assert cost <= tour_cost(dist, greedy(dist, 0, n - 1)) + 1e-9
+
+
+def test_nearest_neighbor_2opt_inf_edge_still_improves():
+    """An unavoidable inf edge must not block finite local improvements.
+
+    Every edge into city 4 is inf, so the total stays inf, but the
+    interior [1, 2] reversal is still an improvement and must be
+    applied (inf-aware prefix sums; naive prefix subtraction would
+    produce nan and reject it).
+    """
+    inf = float('inf')
+    dist = [
+        [0, 1, 2, 50, inf],
+        [1, 0, 1, 5, inf],
+        [2, 1, 0, 10, inf],
+        [50, 5, 10, 0, inf],
+        [inf, inf, inf, inf, 0],
+    ]
+    order, cost = _nearest_neighbor_2opt(dist, 0, 4)
+    assert order == [0, 2, 1, 3, 4]
+    assert cost == inf
+
+
+def test_get_pixel_id_default_dims():
+    """Omitted xdim/ydim fall back to the raster's last two dims."""
+    import xarray as xr
+
+    from xrspatial.pathfinding import _get_pixel_id
+
+    agg = xr.DataArray(np.zeros((5, 5)), dims=['y', 'x'])
+    agg['y'] = np.linspace(4, 0, 5)
+    agg['x'] = np.linspace(0, 4, 5)
+
+    assert _get_pixel_id((4.0, 0.0), agg) == (0, 0)
+    assert _get_pixel_id((0.0, 4.0), agg) == (4, 4)
+    assert (_get_pixel_id((2.0, 3.0), agg)
+            == _get_pixel_id((2.0, 3.0), agg, 'x', 'y'))
+
+
+class TestEmptySurface:
+    """A zero-size surface raises a clear error, not a numpy internal one."""
+
+    @staticmethod
+    def _empty_raster(shape):
+        import xarray as xr
+        r = xr.DataArray(np.ones(shape), dims=('y', 'x'),
+                         attrs={'res': (1.0, 1.0)})
+        r['y'] = np.linspace(shape[0] - 1, 0, shape[0])
+        r['x'] = np.linspace(0, shape[1] - 1, shape[1])
+        return r
+
+    @pytest.mark.parametrize('shape', [(0, 0), (0, 4), (4, 0)])
+    def test_a_star_search_empty_raises(self, shape):
+        r = self._empty_raster(shape)
+        with pytest.raises(ValueError, match='empty'):
+            a_star_search(r, (0.0, 0.0), (1.0, 1.0))
+
+    def test_multi_stop_search_empty_raises(self):
+        r = self._empty_raster((0, 0))
+        with pytest.raises(ValueError, match='empty'):
+            multi_stop_search(r, [(0.0, 0.0), (1.0, 1.0)])
+
+
+@pytest.mark.skipif(not has_dask_array(), reason="Requires dask.Array")
+class TestDaskFrictionFminCache:
+
+    def test_fmin_cached_by_graph_token(self):
+        from xrspatial.pathfinding import _dask_fmin_cache, _dask_friction_fmin
+
+        f = da.from_array(
+            np.array([[2.0, 3.0], [0.5, -1.0]]), chunks=(1, 2))
+        _dask_fmin_cache.clear()
+        assert _dask_friction_fmin(f) == 0.5
+        assert f.name in _dask_fmin_cache
+        # Second call is served from the cache
+        assert _dask_friction_fmin(f) == 0.5
+        assert len(_dask_fmin_cache) == 1
+
+    def test_multi_stop_dask_friction_matches_numpy(self):
+        # Routes 3 segments through the same friction raster, so the
+        # second and third segments hit the f_min cache
+        from xrspatial.pathfinding import _dask_fmin_cache
+
+        data = np.ones((6, 6))
+        friction = np.linspace(1.0, 2.0, 36).reshape(6, 6)
+        agg_np = _make_raster(data)
+        fr_np = _make_raster(friction)
+        waypoints = [(5.0, 0.0), (0.0, 2.0), (5.0, 4.0), (0.0, 5.0)]
+
+        expected = multi_stop_search(agg_np, waypoints, friction=fr_np)
+
+        agg_dask = _make_raster(data)
+        agg_dask.data = da.from_array(agg_dask.data, chunks=(3, 3))
+        fr_dask = _make_raster(friction)
+        fr_dask.data = da.from_array(fr_dask.data, chunks=(3, 3))
+
+        _dask_fmin_cache.clear()
+        result = multi_stop_search(agg_dask, waypoints, friction=fr_dask)
+        # All 3 segments resolved f_min through one cache entry; a
+        # per-call graph token would show up here as extra entries
+        assert len(_dask_fmin_cache) == 1
+
+        np.testing.assert_allclose(
+            result.data.compute(), expected.values, equal_nan=True)
+        assert result.attrs['total_cost'] == pytest.approx(
+            expected.attrs['total_cost'])