ViSR plugin: scan-aware /binned and a prototype /render endpoint

README.md

@@ -26,13 +26,34 @@ tiled server config.
 ### `dls_tiled.visr:visr_router`
 
 `GET /api/v1/binned/{uid}` bins a ViSR scan into a 2-D image on the server and
-returns the mean `RedTotal`, `GreenTotal` and `BlueTotal` in each bin, plus the
-bin edges as `x_limits` and `y_limits`. Positions come from the recorded
-readbacks, or from the ScanSpec setpoints in the run's start document with
-`setpoints=true`. Other query parameters: `x_dim_index`, `y_dim_index`, `width`,
-`height`, `xmin`/`xmax`/`ymin`/`ymax` and the repeatable
+returns the mean `RedTotal`, `GreenTotal` and `BlueTotal` in each bin as
+matrices with one row per y bin and one column per x bin (`null` where no point
+landed), the bin edges as `x_limits` and `y_limits`, and `n_points`/`n_total`,
+the points used so far and the points the scan will have (`null` without a
+spec). Positions come from the recorded readbacks, or from the ScanSpec
+setpoints in the run's start document with `setpoints=true`.
+
+Unless `width`/`height` or a full `xmin`/`xmax`/`ymin`/`ymax` range are given,
+the bins follow the scan's setpoint grid from its spec (one cell per setpoint,
+fixed from the first poll), provided the spec is a grid of at least 3×3 and the
+recorded positions lie on it. Positions that are off the grid by a constant
+amount per axis (as fly scans currently record them) are first shifted onto it
+by the median of each point's position minus its setpoint, and the shift is
+returned as `position_offset`. Otherwise numpy's default 10×10 over the data.
+Other query parameters: `x_dim_index`, `y_dim_index` and the repeatable
 `slice_dim=dim:center:thickness`.
 
+A file still behind the catalog is retried briefly and then answered with 503
+and `Retry-After`. Identical requests share one result for a short while and
+carry an `ETag`, so a repeated poll of an unchanged image gets a `304`.
+
+| Environment variable | Default | |
+|---|---|---|
+| `DLS_TILED_VISR_NOT_READY_RETRIES` | `3` | retries before answering 503 |
+| `DLS_TILED_VISR_NOT_READY_RETRY_DELAY` | `0.2` | seconds between those retries |
+| `DLS_TILED_VISR_CACHE_SECONDS` | `0.5` | how long a result is shared; `0` turns sharing off |
+| `DLS_TILED_VISR_MAX_GRID_CELLS` | `100000` | largest scan grid used by default |
+
 ```yaml
 routers:
   - dls_tiled.visr:visr_router

src/dls_tiled/rendering.py

@@ -0,0 +1,303 @@
+"""Draw scan points into an image for a viewport (prototype).
+
+Each renderer takes the points of one channel (values and positions) and a
+``Viewport`` (a rectangle in data coordinates and an output size in pixels), and
+returns two ``(height, width)`` arrays: the image (NaN where nothing is drawn)
+and the coverage (0 to 1, how much of the pixel is covered by the drawn points).
+Row 0 is the viewport's ``y0`` edge and column 0 its ``x0`` edge.
+
+- ``cells``: the scan's grid cells, each with the mean of the points in it.
+- ``rectangles``: a rectangle of the nominal cell size centred on each recorded
+  position, so the picture follows where the points really were.
+- ``gaussian``: a beam-sized 2-D Gaussian per point, so any gaps the scan left
+  between beam footprints show up as low coverage.
+- ovals (``oval_layout``): a beam-sized oval per point, filled with its value
+  and faintly outlined; overlaps are split between the ovals, not blended.
+"""
+
+from __future__ import annotations
+
+import struct
+import zlib
+from dataclasses import dataclass
+
+import numpy
+
+FWHM_TO_SIGMA = 1 / (2 * numpy.sqrt(2 * numpy.log(2)))
+
+
+@dataclass(frozen=True)
+class Viewport:
+    """A rectangle in data coordinates drawn at ``width`` x ``height`` pixels."""
+
+    x0: float
+    x1: float
+    y0: float
+    y1: float
+    width: int
+    height: int
+
+    @classmethod
+    def fit(cls, x0, x1, y0, y1, size):
+        """The viewport with ``size`` pixels along its longer side, square pixels."""
+        span_x, span_y = abs(x1 - x0) or 1.0, abs(y1 - y0) or 1.0
+        if span_x >= span_y:
+            width, height = size, max(1, round(size * span_y / span_x))
+        else:
+            width, height = max(1, round(size * span_x / span_y)), size
+        return cls(float(x0), float(x1), float(y0), float(y1), int(width), int(height))
+
+    @property
+    def pixel_x(self) -> float:
+        return (self.x1 - self.x0) / self.width
+
+    @property
+    def pixel_y(self) -> float:
+        return (self.y1 - self.y0) / self.height
+
+    def centres(self):
+        """Pixel centres along x (columns) and y (rows)."""
+        xs = self.x0 + (numpy.arange(self.width) + 0.5) * self.pixel_x
+        ys = self.y0 + (numpy.arange(self.height) + 0.5) * self.pixel_y
+        return xs, ys
+
+    def column(self, x):
+        """Fractional column of a data x coordinate."""
+        return (numpy.asarray(x, dtype=float) - self.x0) / self.pixel_x
+
+    def row(self, y):
+        return (numpy.asarray(y, dtype=float) - self.y0) / self.pixel_y
+
+
+def _finish(total, weight, coverage):
+    image = numpy.full(total.shape, numpy.nan)
+    drawn = weight > 0
+    image[drawn] = total[drawn] / weight[drawn]
+    return image, numpy.clip(coverage, 0.0, 1.0)
+
+
+def render_cells(values, x, y, x_edges, y_edges, viewport: Viewport):
+    """The scan's grid cells (mean of the points in each), sampled at each pixel."""
+    values = numpy.asarray(values, dtype=float)
+    counts, _, _ = numpy.histogram2d(x, y, bins=[x_edges, y_edges])
+    sums, _, _ = numpy.histogram2d(x, y, bins=[x_edges, y_edges], weights=values)
+    xs, ys = viewport.centres()
+    col = numpy.searchsorted(x_edges, xs, side="right") - 1
+    row = numpy.searchsorted(y_edges, ys, side="right") - 1
+    inside_x = (col >= 0) & (col < len(x_edges) - 1)
+    inside_y = (row >= 0) & (row < len(y_edges) - 1)
+    total = numpy.zeros((viewport.height, viewport.width))
+    weight = numpy.zeros_like(total)
+    rr, cc = numpy.meshgrid(
+        numpy.where(inside_y)[0], numpy.where(inside_x)[0], indexing="ij"
+    )
+    cell_y, cell_x = row[rr], col[cc]
+    total[rr, cc] = sums[cell_x, cell_y]
+    weight[rr, cc] = counts[cell_x, cell_y]
+    return _finish(total, weight, (weight > 0).astype(float))
+
+
+def render_rectangles(values, x, y, cell_x, cell_y, viewport: Viewport):
+    """A rectangle of ``cell_x`` by ``cell_y`` centred on each point.
+
+    Where rectangles overlap, the pixel takes the mean of their values, weighted
+    by how much of the pixel each one covers.
+    """
+    values = numpy.asarray(values, dtype=float)
+    total = numpy.zeros((viewport.height, viewport.width))
+    weight = numpy.zeros_like(total)
+    for v, px, py in zip(values, viewport.column(x), viewport.row(y), strict=True):
+        if not numpy.isfinite(v):
+            continue
+        cx0, cx1 = (
+            px - cell_x / 2 / viewport.pixel_x,
+            px + cell_x / 2 / viewport.pixel_x,
+        )
+        cy0, cy1 = (
+            py - cell_y / 2 / viewport.pixel_y,
+            py + cell_y / 2 / viewport.pixel_y,
+        )
+        c_lo, c_hi = (
+            max(0, int(numpy.floor(cx0))),
+            min(viewport.width, int(numpy.ceil(cx1))),
+        )
+        r_lo, r_hi = (
+            max(0, int(numpy.floor(cy0))),
+            min(viewport.height, int(numpy.ceil(cy1))),
+        )
+        if c_lo >= c_hi or r_lo >= r_hi:
+            continue
+        cols = numpy.arange(c_lo, c_hi)
+        rows = numpy.arange(r_lo, r_hi)
+        # fraction of each pixel the rectangle covers, along each axis
+        fx = numpy.clip(numpy.minimum(cols + 1, cx1) - numpy.maximum(cols, cx0), 0, 1)
+        fy = numpy.clip(numpy.minimum(rows + 1, cy1) - numpy.maximum(rows, cy0), 0, 1)
+        overlap = numpy.outer(fy, fx)
+        total[r_lo:r_hi, c_lo:c_hi] += overlap * v
+        weight[r_lo:r_hi, c_lo:c_hi] += overlap
+    return _finish(total, weight, weight)
+
+
+def render_gaussian(values, x, y, fwhm_x, fwhm_y, viewport: Viewport, reach=3.0):
+    """A 2-D Gaussian of the beam's size (FWHM in data units) on each point.
+
+    Each pixel takes the mean of the points' values weighted by their Gaussians.
+    Coverage is the summed Gaussian weight, 1 at the centre of an isolated beam,
+    so the space between beam footprints that the scan didn't cover shows up as
+    low coverage rather than being filled in.
+    """
+    values = numpy.asarray(values, dtype=float)
+    sx = fwhm_x * FWHM_TO_SIGMA / viewport.pixel_x  # sigma in pixels
+    sy = fwhm_y * FWHM_TO_SIGMA / viewport.pixel_y
+    total = numpy.zeros((viewport.height, viewport.width))
+    weight = numpy.zeros_like(total)
+    for v, px, py in zip(values, viewport.column(x), viewport.row(y), strict=True):
+        if not numpy.isfinite(v):
+            continue
+        c_lo, c_hi = (
+            max(0, int(px - reach * sx)),
+            min(viewport.width, int(px + reach * sx) + 1),
+        )
+        r_lo, r_hi = (
+            max(0, int(py - reach * sy)),
+            min(viewport.height, int(py + reach * sy) + 1),
+        )
+        if c_lo >= c_hi or r_lo >= r_hi:
+            continue
+        gx = numpy.exp(-0.5 * ((numpy.arange(c_lo, c_hi) + 0.5 - px) / sx) ** 2)
+        gy = numpy.exp(-0.5 * ((numpy.arange(r_lo, r_hi) + 0.5 - py) / sy) ** 2)
+        g = numpy.outer(gy, gx)
+        total[r_lo:r_hi, c_lo:c_hi] += g * v
+        weight[r_lo:r_hi, c_lo:c_hi] += g
+    image, coverage = _finish(total, weight, weight)
+    image[coverage < 0.02] = numpy.nan  # beyond the beams' reach: nothing drawn
+    return image, coverage
+
+
+@dataclass(frozen=True)
+class OvalLayout:
+    """Where each point's oval is drawn in a viewport, independent of any channel.
+
+    ``owner`` is the index of the point each pixel takes its colour from (-1
+    outside every oval) and ``outline`` the opacity (0 to 1) of the white
+    outlines drawn over the colours.
+    """
+
+    owner: numpy.ndarray
+    outline: numpy.ndarray
+
+    def paint(self, values):
+        """One channel's image and coverage: each oval filled with its point's value."""
+        values = numpy.asarray(values, dtype=float)
+        image = numpy.full(self.owner.shape, numpy.nan)
+        inside = self.owner >= 0
+        image[inside] = values[self.owner[inside]]
+        return image, (~numpy.isnan(image)).astype(float)
+
+
+def oval_layout(
+    x, y, width_x, width_y, viewport: Viewport, line_width=1.5, line_opacity=0.35
+) -> OvalLayout:
+    """An oval ``width_x`` by ``width_y`` (data units) centred on each point.
+
+    Where ovals overlap, the overlap is split along the line through the two
+    points where their outlines cross, each side taking the colour of its own
+    oval. For ovals of one size that line is where the two centres are equally
+    far in units of the ovals' radii, so each pixel takes the oval whose
+    centre is nearest in those units; that also splits three-way overlaps
+    consistently.
+
+    Every oval's whole outline is drawn as a faint white line, ``line_width``
+    pixels wide at ``line_opacity`` (0 for none).
+    """
+    a = width_x / 2 / viewport.pixel_x  # semi-axes in pixels
+    b = width_y / 2 / viewport.pixel_y
+    px, py = viewport.column(x), viewport.row(y)
+    shape = (viewport.height, viewport.width)
+    best = numpy.full(shape, numpy.inf)  # the nearest centre, in radii squared
+    owner = numpy.full(shape, -1)
+    outline = numpy.zeros(shape)
+    pad = line_width + 1  # reach past the oval far enough to draw its outline
+    for i in numpy.flatnonzero(numpy.isfinite(px) & numpy.isfinite(py)):
+        c_lo = max(0, int(numpy.floor(px[i] - a - pad)))
+        c_hi = min(viewport.width, int(numpy.ceil(px[i] + a + pad)))
+        r_lo = max(0, int(numpy.floor(py[i] - b - pad)))
+        r_hi = min(viewport.height, int(numpy.ceil(py[i] + b + pad)))
+        if c_lo >= c_hi or r_lo >= r_hi:
+            continue
+        window = (slice(r_lo, r_hi), slice(c_lo, c_hi))
+        dx = (numpy.arange(c_lo, c_hi) + 0.5 - px[i]) / a  # in radii
+        dy = (numpy.arange(r_lo, r_hi) + 0.5 - py[i]) / b
+        q = dy[:, None] ** 2 + dx[None, :] ** 2  # < 1 inside the oval
+        nearer = q < best[window]
+        best[window][nearer] = q[nearer]
+        owner[window][nearer] = i
+        if line_opacity > 0:
+            # distance to the outline in pixels, to first order
+            r = numpy.sqrt(q)
+            with numpy.errstate(divide="ignore", invalid="ignore"):
+                gradient = numpy.hypot(dx[None, :] / a, dy[:, None] / b) / r
+                distance = numpy.where(r > 0, (r - 1) / gradient, numpy.inf)
+            line = numpy.clip(line_width / 2 + 0.5 - numpy.abs(distance), 0.0, 1.0)
+            outline[window] = numpy.maximum(outline[window], line * line_opacity)
+    return OvalLayout(numpy.where(best < 1, owner, -1), outline)
+
+
+def to_rows(image) -> list[list[float | None]]:
+    """An image as JSON rows, ``None`` where nothing was drawn or the value isn't
+    a number JSON can carry (NaN or infinite)."""
+    out = image.astype(object)
+    out[~numpy.isfinite(image)] = None
+    return out.tolist()
+
+
+CHANNEL_COLOUR = {"RedTotal": 0, "GreenTotal": 1, "BlueTotal": 2}
+
+
+def to_png(image, coverage, colour=0, vmin=None, vmax=None, outline=None) -> bytes:
+    """An RGBA PNG of one channel, tinted (0/1/2 = red/green/blue), transparent
+    where nothing was drawn, with an optional white ``outline`` layer (opacity
+    0 to 1) drawn over it. Row 0 of the image is the top row of the PNG."""
+    valid = image[~numpy.isnan(image)]
+    lo = float(valid.min()) if vmin is None and valid.size else (vmin or 0.0)
+    hi = float(valid.max()) if vmax is None and valid.size else (vmax or 1.0)
+    scaled = numpy.zeros(image.shape)
+    if hi > lo:
+        scaled = numpy.clip((numpy.nan_to_num(image, nan=lo) - lo) / (hi - lo), 0, 1)
+    rgb = numpy.zeros(image.shape + (3,))
+    rgb[..., colour] = scaled
+    alpha = numpy.where(numpy.isnan(image), 0.0, coverage)
+    if outline is not None:
+        # white over the colour, where the colour itself may be transparent
+        total = outline + alpha * (1 - outline)
+        with numpy.errstate(divide="ignore", invalid="ignore"):
+            rgb = numpy.where(
+                total[..., None] > 0,
+                (outline[..., None] + rgb * (alpha * (1 - outline))[..., None])
+                / total[..., None],
+                0.0,
+            )
+        alpha = total
+    rgba = numpy.zeros(image.shape + (4,), dtype=numpy.uint8)
+    rgba[..., :3] = numpy.clip(numpy.round(rgb * 255), 0, 255)
+    rgba[..., 3] = numpy.clip(numpy.round(alpha * 255), 0, 255)
+    return _encode_png(rgba)
+
+
+def _encode_png(rgba) -> bytes:
+    height, width, _ = rgba.shape
+    raw = b"".join(b"\x00" + rgba[r].tobytes() for r in range(height))
+
+    def chunk(kind: bytes, payload: bytes) -> bytes:
+        body = kind + payload
+        return (
+            struct.pack(">I", len(payload)) + body + struct.pack(">I", zlib.crc32(body))
+        )
+
+    header = struct.pack(">IIBBBBB", width, height, 8, 6, 0, 0, 0)
+    return (
+        b"\x89PNG\r\n\x1a\n"
+        + chunk(b"IHDR", header)
+        + chunk(b"IDAT", zlib.compress(raw, 6))
+        + chunk(b"IEND", b"")
+    )