Address review: use scikit-image, vectorize TileGrid, remove double ensure_f32_2d

- Replace numba @njit convolutions with skimage.filters.sobel_h/sobel_v
  (tenengrad), skimage.filters.laplace (laplacian variance), and np.var
  (population variance). Removes numba dependency from sharpness metrics.
- Vectorize TileGrid: replace itertools.product loops with numpy
  broadcasting for indices/bounds and shapely.box for batch polygon
  creation.
- Remove redundant _ensure_f32_2d calls from metric registry lambdas
  (each metric function already handles its own input validation).

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

Author: timtreis

src/squidpy/experimental/im/_qc_metrics.py

@@ -7,7 +7,6 @@
 import numpy as np
 
 from squidpy.experimental.im._sharpness_metrics import (
-    _ensure_f32_2d,
     _fft_high_freq_energy,
     _haar_wavelet_energy,
     _laplacian_variance,
@@ -186,9 +185,9 @@ def _tissue_fraction(block: np.ndarray) -> np.ndarray:
 
 _METRIC_REGISTRY: dict[QCMetric, tuple[InputKind, MetricFn]] = {
     # Sharpness (grayscale)
-    QCMetric.TENENGRAD: (InputKind.GRAYSCALE, lambda a: _tenengrad_mean(_ensure_f32_2d(a))),
-    QCMetric.VAR_OF_LAPLACIAN: (InputKind.GRAYSCALE, lambda a: _laplacian_variance(_ensure_f32_2d(a))),
-    QCMetric.VARIANCE: (InputKind.GRAYSCALE, lambda a: _pop_variance(_ensure_f32_2d(a))),
+    QCMetric.TENENGRAD: (InputKind.GRAYSCALE, _tenengrad_mean),
+    QCMetric.VAR_OF_LAPLACIAN: (InputKind.GRAYSCALE, _laplacian_variance),
+    QCMetric.VARIANCE: (InputKind.GRAYSCALE, _pop_variance),
     QCMetric.FFT_HIGH_FREQ_ENERGY: (InputKind.GRAYSCALE, _fft_high_freq_energy),
     QCMetric.HAAR_WAVELET_ENERGY: (InputKind.GRAYSCALE, _haar_wavelet_energy),
     # Intensity (grayscale)

src/squidpy/experimental/im/_sharpness_metrics.py

@@ -1,8 +1,8 @@
 from __future__ import annotations
 
 import numpy as np
-from numba import njit
 from scipy.fft import fft2, fftfreq
+from skimage.filters import laplace, sobel_h, sobel_v
 
 
 def _ensure_f32_2d(x: np.ndarray) -> np.ndarray:
@@ -11,79 +11,25 @@ def _ensure_f32_2d(x: np.ndarray) -> np.ndarray:
     return np.ascontiguousarray(x.astype(np.float32, copy=False))
 
 
-@njit(cache=True, fastmath=True)
-def _clamp(v: int, lo: int, hi: int) -> int:
-    if v < lo:
-        return lo
-    if v > hi:
-        return hi
-    return v
+def _tenengrad_mean(block: np.ndarray) -> np.ndarray:
+    """Mean Tenengrad energy (sum of squared Sobel gradients)."""
+    b = _ensure_f32_2d(block)
+    energy = sobel_h(b) ** 2 + sobel_v(b) ** 2
+    return np.array([[float(energy.mean())]], dtype=np.float32)
 
 
-@njit(cache=True, fastmath=True)
-def _tenengrad_mean(block: np.ndarray) -> np.ndarray:
-    """Mean Tenengrad energy using Sobel 3x3."""
-    h, w = block.shape
-    gxk = np.array([[-1.0, 0.0, 1.0], [-2.0, 0.0, 2.0], [-1.0, 0.0, 1.0]], dtype=np.float32)
-    gyk = np.array([[1.0, 2.0, 1.0], [0.0, 0.0, 0.0], [-1.0, -2.0, -1.0]], dtype=np.float32)
-    s = 0.0
-    for i in range(h):
-        for j in range(w):
-            gx = 0.0
-            gy = 0.0
-            for di in range(-1, 2):
-                for dj in range(-1, 2):
-                    ii = _clamp(i + di, 0, h - 1)
-                    jj = _clamp(j + dj, 0, w - 1)
-                    v = block[ii, jj]
-                    gx += gxk[di + 1, dj + 1] * v
-                    gy += gyk[di + 1, dj + 1] * v
-            s += gx * gx + gy * gy
-    mean_val = s / (h * w)
-    return np.array([[mean_val]], dtype=np.float32)
-
-
-@njit(cache=True, fastmath=True)
 def _laplacian_variance(block: np.ndarray) -> np.ndarray:
     """Population variance of Laplacian response."""
-    h, w = block.shape
-    lk = np.array([[0.0, 1.0, 0.0], [1.0, -4.0, 1.0], [0.0, 1.0, 0.0]], dtype=np.float32)
-    n = h * w
-    s = 0.0
-    s2 = 0.0
-    for i in range(h):
-        for j in range(w):
-            y = 0.0
-            for di in range(-1, 2):
-                for dj in range(-1, 2):
-                    ii = _clamp(i + di, 0, h - 1)
-                    jj = _clamp(j + dj, 0, w - 1)
-                    y += lk[di + 1, dj + 1] * block[ii, jj]
-            s += y
-            s2 += y * y
-    mean = s / n
-    # var = E[y^2] - (E[y])^2
-    var = (s2 / n) - (mean * mean)
-    var_val = var if var > 0.0 else 0.0
-    return np.array([[var_val]], dtype=np.float32)
-
-
-@njit(cache=True, fastmath=True)
+    b = _ensure_f32_2d(block)
+    lap = laplace(b)
+    var_val = float(np.var(lap))
+    return np.array([[max(var_val, 0.0)]], dtype=np.float32)
+
+
 def _pop_variance(block: np.ndarray) -> np.ndarray:
     """Population variance of pixel intensities."""
-    h, w = block.shape
-    n = h * w
-    s = 0.0
-    s2 = 0.0
-    for i in range(h):
-        for j in range(w):
-            v = block[i, j]
-            s += v
-            s2 += v * v
-    mean = s / n
-    var = (s2 / n) - (mean * mean)
-    var_val = var if var > 0.0 else 0.0
-    return np.array([[var_val]], dtype=np.float32)
+    b = _ensure_f32_2d(block)
+    return np.array([[float(np.var(b))]], dtype=np.float32)
 
 
 def _fft_high_freq_energy(block: np.ndarray) -> np.ndarray:

src/squidpy/experimental/im/_utils.py

@@ -1,13 +1,12 @@
 from __future__ import annotations
 
-import itertools
 from typing import TYPE_CHECKING, Any, Literal
 
 import dask.array as da
 import geopandas as gpd
 import numpy as np
 import xarray as xr
-from shapely.geometry import Polygon
+from shapely import box
 from spatialdata._logging import logger
 from spatialdata.models import ShapesModel
 from spatialdata.transformations import get_transformation, set_transformation
@@ -46,13 +45,13 @@ def __init__(
         total_w_needed = self.W - grid_start_x
         self.tiles_y = (total_h_needed + self.ty - 1) // self.ty
         self.tiles_x = (total_w_needed + self.tx - 1) // self.tx
-        # Cache immutable derived values
-        self._indices = np.array(
-            [[iy, ix] for iy, ix in itertools.product(range(self.tiles_y), range(self.tiles_x))], dtype=int
-        )
-        self._names = [f"tile_x{ix}_y{iy}" for iy, ix in itertools.product(range(self.tiles_y), range(self.tiles_x))]
-        self._bounds = self._compute_bounds()
-        self._centroids_polys = self._compute_centroids_and_polygons()
+        # Cache immutable derived values (vectorized)
+        iy = np.repeat(np.arange(self.tiles_y), self.tiles_x)
+        ix = np.tile(np.arange(self.tiles_x), self.tiles_y)
+        self._indices = np.column_stack([iy, ix])
+        self._names = [f"tile_x{x}_y{y}" for y, x in zip(iy, ix, strict=True)]
+        self._bounds = self._compute_bounds(iy, ix)
+        self._centroids, self._polys = self._compute_centroids_and_polygons()
 
     def indices(self) -> np.ndarray:
         return self._indices
@@ -63,33 +62,29 @@ def names(self) -> list[str]:
     def bounds(self) -> np.ndarray:
         return self._bounds
 
-    def _compute_bounds(self) -> np.ndarray:
-        b: list[list[int]] = []
-        for iy, ix in itertools.product(range(self.tiles_y), range(self.tiles_x)):
-            y0 = iy * self.ty + self.offset_y
-            x0 = ix * self.tx + self.offset_x
-            y1 = ((iy + 1) * self.ty + self.offset_y) if iy < self.tiles_y - 1 else self.H
-            x1 = ((ix + 1) * self.tx + self.offset_x) if ix < self.tiles_x - 1 else self.W
-            # Clamp bounds to image dimensions
-            y0 = max(0, min(y0, self.H))
-            x0 = max(0, min(x0, self.W))
-            y1 = max(0, min(y1, self.H))
-            x1 = max(0, min(x1, self.W))
-            b.append([y0, x0, y1, x1])
-        return np.array(b, dtype=int)
-
-    def centroids_and_polygons(self) -> tuple[np.ndarray, list[Polygon]]:
-        return self._centroids_polys
-
-    def _compute_centroids_and_polygons(self) -> tuple[np.ndarray, list[Polygon]]:
-        cents: list[list[float]] = []
-        polys: list[Polygon] = []
-        for y0, x0, y1, x1 in self._bounds:
-            cy = (y0 + y1) / 2
-            cx = (x0 + x1) / 2
-            cents.append([cy, cx])
-            polys.append(Polygon([(x0, y0), (x1, y0), (x1, y1), (x0, y1), (x0, y0)]))
-        return np.array(cents, dtype=float), polys
+    def _compute_bounds(self, iy: np.ndarray, ix: np.ndarray) -> np.ndarray:
+        y0 = iy * self.ty + self.offset_y
+        x0 = ix * self.tx + self.offset_x
+        y1 = (iy + 1) * self.ty + self.offset_y
+        x1 = (ix + 1) * self.tx + self.offset_x
+        # Last row/column extends to image edge
+        y1[iy == self.tiles_y - 1] = self.H
+        x1[ix == self.tiles_x - 1] = self.W
+        # Clamp to image dimensions
+        y0 = np.clip(y0, 0, self.H)
+        x0 = np.clip(x0, 0, self.W)
+        y1 = np.clip(y1, 0, self.H)
+        x1 = np.clip(x1, 0, self.W)
+        return np.column_stack([y0, x0, y1, x1]).astype(int)
+
+    def centroids_and_polygons(self) -> tuple[np.ndarray, list]:
+        return self._centroids, self._polys
+
+    def _compute_centroids_and_polygons(self) -> tuple[np.ndarray, list]:
+        y0, x0, y1, x1 = self._bounds[:, 0], self._bounds[:, 1], self._bounds[:, 2], self._bounds[:, 3]
+        centroids = np.column_stack([(y0 + y1) / 2.0, (x0 + x1) / 2.0])
+        polys = list(box(x0, y0, x1, y1))
+        return centroids, polys
 
     def rechunk_and_pad(self, arr_yx: da.Array) -> da.Array:
         if arr_yx.ndim != 2: