@@ -167,13 +167,41 @@ def __init__(self, mask: Mask2D, sub_size: Union[int, Array2D]):
         ):
             self.segment_ids[start:end] = seg_id
 
-    @property
+    @cached_property
     def sub_is_uniform(self) -> bool:
         """
         Returns True if the sub_size is uniform across all pixels in the mask.
+
+        The `sub_size` is fixed when the `OverSampler` is created, so this is cached: it is evaluated once per
+        instance instead of on every call to `binned_array_2d_from` (which is called once per light profile
+        evaluation, and therefore many hundreds of times per likelihood evaluation).
         """
         return np.all(np.isclose(self.sub_size, self.sub_size[0]))
 
+    @cached_property
+    def binned_counts(self) -> np.ndarray:
+        """
+        The number of sub-pixels which are binned up into every pixel of the mask, used as the divisor when
+        binning an over sampled array via `binned_array_2d_from` for a non-uniform `sub_size`.
+
+        Segments which contain no sub-pixels have their count set to 1, so that the division by this array does
+        not raise a divide-by-zero (the corresponding sum is zero, so the binned value is zero either way).
+
+        This is a function of the `segment_ids` alone, which are fixed when the `OverSampler` is created, so it is
+        computed once per instance instead of on every call to `binned_array_2d_from`.
+
+        The returned array is read-only, because it is shared by every call and must never be mutated in place.
+        """
+        counts = np.bincount(self.segment_ids, minlength=self.mask.pixels_in_mask)
+
+        # Avoid division by zero (on a copy, because `np.bincount` output is not reused elsewhere but the cached
+        # array below is handed out to every caller).
+
+        counts = np.maximum(counts, 1)
+        counts.flags.writeable = False
+
+        return counts
+
     def tree_flatten(self):
         return (self.mask, self.sub_size), ()
 
@@ -263,13 +291,8 @@ def binned_array_2d_from(self, array: Array2D, xp=np) -> "Array2D":
                     self.segment_ids, weights=array, minlength=self.mask.pixels_in_mask
                 )
 
-                # Count number of items per segment
-                counts = np.bincount(
-                    self.segment_ids, minlength=self.mask.pixels_in_mask
-                )
-
-                # Avoid division by zero
-                counts[counts == 0] = 1
+                # Count number of items per segment (cached, as it depends only on the segment ids)
+                counts = self.binned_counts
 
             binned_array_2d = sums / counts
 
 
@@ -107,3 +107,94 @@ def test__slim_index_for_sub_slim_index():
     )
 
     assert (over_sampling.slim_for_sub_slim == slim_index_for_sub_slim_index_util).all()
+
+
+def test__over_sampler_binned_array_2d__non_uniform_counts_are_cached_and_not_mutated():
+    """
+    The divisor used to bin a non-uniform over sampled array depends only on the `segment_ids`, so it is cached on
+    the `OverSampler`. Repeated calls must therefore be bit-identical and must never mutate the cached array.
+    """
+    mask = aa.Mask2D(
+        mask=[[False, False], [True, True]],
+        pixel_scales=1.0,
+    )
+
+    over_sampler = aa.OverSampler(
+        mask=mask, sub_size=aa.Array2D(values=[1, 2], mask=mask)
+    )
+
+    assert not over_sampler.sub_is_uniform
+
+    arr = np.array([1.0, 5.0, 7.0, 10.0, 10.0])
+
+    binned_0 = over_sampler.binned_array_2d_from(array=arr)
+
+    counts = over_sampler.binned_counts
+
+    assert counts == pytest.approx(np.array([1, 4]), 1.0e-4)
+    assert not counts.flags.writeable
+
+    binned_1 = over_sampler.binned_array_2d_from(array=arr)
+    binned_2 = over_sampler.binned_array_2d_from(array=2.0 * arr)
+
+    assert np.array_equal(np.array(binned_0.slim), np.array(binned_1.slim))
+    assert np.array_equal(np.array(binned_2.slim), 2.0 * np.array(binned_0.slim))
+
+    # The cached array is handed to every caller, so it must be unchanged after the calls above.
+
+    assert over_sampler.binned_counts is counts
+    assert np.array_equal(counts, np.array([1, 4]))
+
+
+def test__over_sampler_binned_array_2d__zero_count_segments_do_not_divide_by_zero():
+    """
+    A `sub_size` of zero produces a segment with no sub-pixels, whose count is guarded to 1 so the binned value is
+    zero rather than a divide-by-zero NaN.
+    """
+    mask = aa.Mask2D(
+        mask=[[False, False], [True, True]],
+        pixel_scales=1.0,
+    )
+
+    with np.errstate(divide="ignore"):
+        over_sampler = aa.OverSampler(
+            mask=mask, sub_size=aa.Array2D(values=[0, 2], mask=mask)
+        )
+
+    arr = np.array([1.0, 2.0, 3.0, 4.0])
+
+    binned = over_sampler.binned_array_2d_from(array=arr)
+
+    assert over_sampler.binned_counts == pytest.approx(np.array([1, 4]), 1.0e-4)
+    assert binned.slim == pytest.approx(np.array([0.0, 2.5]), 1.0e-4)
+
+
+def test__over_sampler__pickles_with_cached_properties():
+    """
+    `OverSampler` instances are pickled to Nautilus worker processes, so the cached divisor must survive a
+    round trip.
+    """
+    import pickle
+
+    mask = aa.Mask2D(
+        mask=[[False, False], [True, True]],
+        pixel_scales=1.0,
+    )
+
+    over_sampler = aa.OverSampler(
+        mask=mask, sub_size=aa.Array2D(values=[1, 2], mask=mask)
+    )
+
+    arr = np.array([1.0, 5.0, 7.0, 10.0, 10.0])
+
+    binned = over_sampler.binned_array_2d_from(array=arr)
+
+    over_sampler_pickled = pickle.loads(pickle.dumps(over_sampler))
+
+    assert np.array_equal(
+        over_sampler_pickled.binned_counts, over_sampler.binned_counts
+    )
+    assert np.array_equal(
+        np.array(over_sampler_pickled.binned_array_2d_from(array=arr).slim),
+        np.array(binned.slim),
+    )