Add pydeserver --twin: serve frames from the de-twin digital twin

deapi/simulated_server/initialize_server.py

@@ -28,13 +28,42 @@ def _recv_exact(conn, n):
     return buf
 
 
+def serve_twin(port, twin_args=()):
+    """Serve frames rendered by the de-twin digital twin instead of the built-in fake data.
+
+    Needs the optional ``twin`` extra (``pip install "deapi[twin]"``). ``twin_args`` are
+    passed to the twin's server, e.g. ``["--specimen", "Apoferritin in ice", "--camera",
+    "Celeritas", "--soap-port", "5002"]``; see ``python -m de_twin.faces.deapi_server --help``.
+    """
+    try:
+        from de_twin.faces import deapi_server
+    except ImportError:
+        sys.stderr.write(
+            'pydeserver --twin needs the digital twin: pip install "deapi[twin]"\n'
+        )
+        return 2
+    return deapi_server.main([str(port), *twin_args])
+
+
 # Defining main function
 def main(port=13240):
     parser = argparse.ArgumentParser()
     parser.add_argument("--port", type=int, help="Port to listen on")
-    args, _ = parser.parse_known_args()
+    parser.add_argument(
+        "--twin",
+        action="store_true",
+        help='serve frames from the de-twin digital twin (pip install "deapi[twin]"); '
+        "other options are passed to the twin, e.g. --specimen, --camera, --soap-port",
+    )
+    args, rest = parser.parse_known_args()
     if args.port:
         port = args.port
+    if args.twin:
+        # the twin may call back into this loop (its --deapi-loop option): don't recurse
+        sys.argv = [a for a in sys.argv if a != "--twin"]
+        if rest and rest[0].isdigit():  # the positional port, already in ``port``
+            rest = rest[1:]
+        return serve_twin(port, rest)
 
     HOST = "127.0.0.1"  # Standard loopback interface address (localhost)
     PORT = port  # Port to listen on (non-privileged ports are > 1023)

deapi/tests/test_fake_server/test_twin_server.py

@@ -0,0 +1,96 @@
+"""``pydeserver --twin``: the simulated server backed by the de-twin digital twin."""
+
+import socket
+import subprocess
+import sys
+import threading
+import time
+
+import numpy as np
+import pytest
+
+from deapi.client import Client
+from deapi.simulated_server import initialize_server
+
+
+def _free_port():
+    with socket.socket() as s:
+        s.bind(("127.0.0.1", 0))
+        return s.getsockname()[1]
+
+
+def test_twin_flag_without_the_twin_installed_explains_how_to_get_it(
+    monkeypatch, capsys
+):
+    monkeypatch.setitem(sys.modules, "de_twin", None)  # import de_twin -> ImportError
+    monkeypatch.setitem(sys.modules, "de_twin.faces", None)
+    monkeypatch.setattr(sys, "argv", ["pydeserver", "--twin"])
+    assert initialize_server.main(_free_port()) == 2
+    assert 'pip install "deapi[twin]"' in capsys.readouterr().err
+
+
+def test_twin_server_serves_twin_frames():
+    pytest.importorskip("de_twin")
+    port = _free_port()
+    proc = subprocess.Popen(
+        [
+            sys.executable,
+            "-u",
+            "-m",
+            "deapi.simulated_server.initialize_server",
+            str(port),
+            "--twin",
+            "--camera",
+            "DESim",
+            "--specimen",
+            "Dense Au on holey C",
+        ],
+        stdout=subprocess.PIPE,
+        stderr=subprocess.STDOUT,
+        text=True,
+    )
+    lines = []
+    started = threading.Event()
+
+    def read():
+        for line in proc.stdout:
+            lines.append(line)
+            if "started" in line:
+                started.set()
+
+    threading.Thread(target=read, daemon=True).start()
+    try:
+        assert started.wait(120), "".join(lines)
+        client = Client()
+        client.usingMmf = False
+        client.connect(port=port)
+        assert client["Sensor Size X (pixels)"] == 1024  # the twin's DESim model
+        client["Frame Count"] = 2
+        client.start_acquisition(1)
+        deadline = time.time() + 60
+        while client.acquiring and time.time() < deadline:
+            time.sleep(0.05)
+        image, *_ = client.get_result("singleframe_integrated")
+        assert image.shape == (1024, 1024)
+        assert np.asarray(image).std() > 0  # a rendered specimen, not a constant
+        assert client["Instrument Project Magnification"]  # the twin's column metadata
+        client.disconnect()
+    finally:
+        proc.terminate()
+        proc.wait(10)
+
+
+def test_twin_options_are_passed_through(monkeypatch):
+    seen = {}
+    monkeypatch.setattr(
+        initialize_server,
+        "serve_twin",
+        lambda port, args: seen.update(port=port, args=args),
+    )
+    monkeypatch.setattr(
+        sys,
+        "argv",
+        ["pydeserver", "13241", "--twin", "--soap-port", "5002", "--camera", "DE16"],
+    )
+    initialize_server.main(13241)
+    assert seen == {"port": 13241, "args": ["--soap-port", "5002", "--camera", "DE16"]}

doc/help/pyDEServer.rst

@@ -62,4 +62,27 @@ The pyDEServer "Cheats" in a couple of ways:
    points of interest.
 
 This is implemented with the `BaseFakeData` class. Which implements a `__getitem__` method that returns the
-data at the index in the navigation data. This can be used to return a single frame or a set of frames.
+data at the index in the navigation data. This can be used to return a single frame or a set of frames.
+Realistic data from the digital twin
+------------------------------------
+
+For data that behaves like a real instrument, the pyDEServer can serve frames rendered by
+`de-twin <https://github.com/directelectron/de-twin>`_, a digital twin of a Direct Electron
+camera on a TEM (column, specimen, in-situ holder and detector; Python 3.10+). Install the optional extra and
+start the server with ``--twin``:
+
+.. code-block::
+
+    pip install "deapi[twin]"
+    pydeserver --port 13241 --twin --camera DE16 --specimen "Dense Au on holey C"
+
+Clients connect exactly as before (``client.usingMmf = False``). The twin supplies:
+
+* raw detector frames with dark offset, noise and gain structure, and references that
+  work like DE-Server's;
+* images that follow the simulated column (stage, magnification, defocus);
+* the ``Instrument ...`` metadata properties.
+
+Other options are passed to the twin. For example, ``--soap-port 5002`` also serves the
+twin's microscope as a DE-TEM-Channel, and ``--seed``, ``--holder`` and ``--time-scale`` are
+available too; see ``python -m de_twin.faces.deapi_server --help``.

examples/digital_twin/README.rst

@@ -0,0 +1,6 @@
+Digital twin
+------------
+
+Examples that run against the `de-twin <https://github.com/directelectron/de-twin>`_ digital
+twin (``pydeserver --twin``): a simulated camera, microscope and specimen that respond to each
+other like a real instrument.

examples/digital_twin/imaging_holes_with_the_digital_twin_sgskip.py

@@ -0,0 +1,147 @@
+"""
+Imaging holes on a grid with the digital twin
+=============================================
+
+A small automation loop that drives the microscope and the camera together, the way it
+would run on an instrument:
+
+1. take a low-magnification atlas of a grid square,
+2. find the holes in the holey carbon film,
+3. move the stage to each hole with ``de_microscope`` (the DE-TEM-Channel client) and take
+   an image at higher magnification.
+
+Here the microscope and the camera are the `de-twin <https://github.com/directelectron/de-twin>`_
+digital twin, so the example runs on any computer. Start it in a terminal first; it serves
+a DE Server (deapi) on port 13250 and a DE-TEM-Channel on port 5002 for the same simulated
+instrument, so stage moves change the images:
+
+.. code-block:: bash
+
+    pip install "deapi[twin]"
+    pydeserver --port 13250 --twin --soap-port 5002 --camera DE16 --specimen "Dense Au on holey C"
+
+Against a real instrument only the addresses change: DE-Server's port and the
+DE-TEM-Channel host.
+
+.. image:: /_static/digital_twin_holes.png
+   :alt: Atlas with the holes found, and an image of each hole
+"""
+
+import time
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy import ndimage
+
+import deapi
+from de_microscope import Microscope
+
+client = deapi.Client()
+client.usingMmf = False  # the twin sends images over the socket
+client.connect(port=13250)
+scope = Microscope(host="127.0.0.1", port=5002)
+
+
+def acquire(frames=10, fps=40):
+    """One integrated, dark/gain-corrected image (like DE-MC's 'Single')."""
+    client["Frames Per Second"] = fps
+    client["Frame Count"] = frames
+    client.start_acquisition(1)
+    while client.acquiring:
+        time.sleep(0.05)
+    image, *_ = client.get_result("singleframe_integrated")
+    return np.asarray(image, dtype=float)
+
+
+def pixel_size_um():
+    """Specimen pixel size reported by the server for the current magnification."""
+    return client["Specimen Pixel Size X (nanometers)"] / 1000.0
+
+
+# %%
+# A low-magnification atlas
+# -------------------------
+# Spread the beam and go to low magnification so a few holes of the holey carbon fit in the
+# field of view. Holes are where the beam passes through vacuum, so they are the brightest
+# regions of the image.
+
+scope["Intensity"] = 0.8
+scope.set("Magnification", 2000, wait=True)
+scope.set(
+    "StagePosition", {"x": 40.0, "y": 0.0}, wait=True
+)  # into a grid square, off the bar
+start = scope["StagePosition"]
+atlas = acquire()
+atlas_px_um = pixel_size_um()
+
+# %%
+# Find the holes
+# --------------
+# Smooth, threshold the bright regions, and keep blobs of a plausible hole size
+# (about 2 um across on this film) that lie wholly inside the atlas; a hole cut by the
+# image edge would give a biased centre. Each hole's offset from the image centre,
+# converted to micrometres, is how far the stage has to move to centre it.
+
+smooth = ndimage.gaussian_filter(atlas, 8)
+bright = smooth > np.percentile(smooth, 80)
+labels, n = ndimage.label(bright)
+areas = ndimage.sum(np.ones_like(atlas), labels, range(1, n + 1)) * atlas_px_um**2
+centres = ndimage.center_of_mass(bright, labels, range(1, n + 1))
+ny, nx = atlas.shape
+margin = 1.5 / atlas_px_um  # pixels: a hole radius plus a little
+holes = [
+    ((c - nx / 2) * atlas_px_um, (r - ny / 2) * atlas_px_um, (r, c))
+    for (r, c), a in zip(centres, areas)
+    if 1.0 < a < 8.0  # um^2: about one hole
+    and margin < r < ny - margin
+    and margin < c < nx - margin
+]
+print(f"found {len(holes)} holes")
+
+# %%
+# Visit each hole
+# ---------------
+# Stage moves carry the specimen with them, so a feature at +dx in the image is centred by
+# moving the stage by -dx. On a real column, calibrate the sign and rotation once
+# (for example with ``de_microscope``'s stage calibration) and apply them here.
+
+scope.set("Magnification", 8000, wait=True)
+scope["Intensity"] = 0.6  # converge the beam again for the close-ups
+hole_images = []
+for dx_um, dy_um, _ in holes[:4]:
+    scope.set(
+        "StagePosition", {"x": start["x"] - dx_um, "y": start["y"] - dy_um}, wait=True
+    )
+    hole_images.append(acquire(frames=40))  # 1 s
+scope.set("StagePosition", {"x": start["x"], "y": start["y"]}, wait=True)
+
+# %%
+# Plot the atlas and the holes
+# ----------------------------
+
+
+def show(ax, img):
+    # the carbon film only dims the beam by ~25 %, while the gold particles are nearly
+    # black: set the grey range from the film and the holes, not the particles
+    lo, hi = np.percentile(img, [10, 99.8])
+    ax.imshow(img[::4, ::4], cmap="gray", vmin=lo, vmax=hi)
+
+
+fig, axes = plt.subplots(
+    1, 1 + len(hole_images), figsize=(3.2 * (1 + len(hole_images)), 3.4)
+)
+axes = np.atleast_1d(axes)
+show(axes[0], atlas)
+for i, (_, _, (r, c)) in enumerate(holes[: len(hole_images)]):
+    axes[0].plot(c / 4, r / 4, "o", mfc="none", mec="C1", ms=14)
+    axes[0].text(c / 4, r / 4, str(i + 1), color="C1", ha="center", va="center")
+axes[0].set_title("atlas, 2000x")
+for i, (ax, img) in enumerate(zip(axes[1:], hole_images)):
+    show(ax, img)
+    ax.set_title(f"hole {i + 1}, 8000x")
+for ax in axes:
+    ax.axis("off")
+fig.tight_layout()
+plt.show()
+
+client.disconnect()

pyproject.toml

@@ -60,6 +60,10 @@ tests = [
 service = [
     "pyqt6",
 ]
+# Serve frames from the de-twin digital twin: `pydeserver --twin`
+twin = [
+    "de-twin>=0.1.1; python_version >= '3.10'",
+]
 doc = [
     "sphinx",
     "pydata_sphinx_theme",

upcoming_changes/59.doc.rst

@@ -0,0 +1 @@
+Added an example that moves the stage with ``de_microscope`` and images holes in a grid, run against the digital twin.

upcoming_changes/59.new_feature.rst

@@ -0,0 +1 @@
+``pydeserver --twin`` serves frames rendered by the `de-twin <https://github.com/directelectron/de-twin>`_ digital twin instead of the built-in fake data: a simulated microscope, specimen and detector, with a DE-TEM-Channel (``--soap-port``) so stage and optics changes show up in the images. Install it with ``pip install "deapi[twin]"``.