Surface detection and faster z-stacks

acq4/devices/Camera/CameraInterface.py

@@ -346,8 +346,10 @@ def boundingRect(self):
         return self.cam.getBoundary().boundingRect()
 
     def takeImage(self, closeShutter=None):
+        # TODO this is maybe unused
         # closeShutter is used for laser scanning devices; we can ignore it here.
-        return self.getDevice().acquireFrames(1, stack=False)
+        with self.getDevice().run():
+            return self.getDevice().acquireFrames(1).getResult()[0]
 
 
 class CameraItemGroup(DeviceTreeItemGroup):

acq4/devices/Camera/taskGUI.py

@@ -85,7 +85,7 @@ def generateTask(self, params=None):
             task['pushState'] = None
             task['popState'] = None
         if state['fixedFrameEnabled']:
-            task['minFrames'] = state['minFrames']
+            task['minFrames'] = int(state['minFrames'])
         return task
 
     def taskSequenceStarted(self):

acq4/devices/DAQGeneric/DaqChannelGui.py

@@ -287,6 +287,8 @@ def updateHolding(self):
         hv = self.getHoldingValue()
         if hv is not None:
             if not self.ui.holdingCheck.isChecked():
+                if self.config['type'] == 'do':
+                    hv = int(hv)
                 self.ui.holdingSpin.setValue(hv)
             self.ui.waveGeneratorWidget.setOffset(hv)
 

acq4/devices/Device.py

@@ -3,13 +3,15 @@
 
 import os
 import traceback
-import weakref
+from contextlib import contextmanager
+from typing import Optional
 
 import acq4
 from acq4.Interfaces import InterfaceMixin
 from acq4.util import Qt
 from acq4.util.Mutex import Mutex
 from acq4.util.debug import printExc
+from acq4.util.optional_weakref import Weakref
 
 
 class Device(InterfaceMixin, Qt.QObject):  # QObject calls super, which is disastrous if not last in the MRO
@@ -30,7 +32,7 @@ def __init__(self, deviceManager: acq4.Manager.Manager, config: dict, name: str)
         self._lock_tb_ = None
         self.dm = deviceManager
         self.dm.declareInterface(name, ['device'], self)
-        Device._deviceCreationOrder.append(weakref.ref(self))
+        Device._deviceCreationOrder.append(Weakref(self))
         self._name = name
 
     def name(self):
@@ -86,6 +88,14 @@ def appendConfigFile(self, data, filename):
         fileName = os.path.join(self.configPath(), filename)
         return self.dm.appendConfigFile(data, fileName)
 
+    @contextmanager
+    def reserved(self):
+        self.reserve()
+        try:
+            yield
+        finally:
+            self.release()
+
     def reserve(self, block=True, timeout=20):
         """Reserve this device globally.
 
@@ -160,7 +170,7 @@ def __init__(self, dev, cmd, parentTask):
         operating synchronously.
         """
         self.dev = dev
-        self.__parentTask = weakref.ref(parentTask)
+        self.__parentTask = Weakref(parentTask)
 
     def parentTask(self):
         return self.__parentTask()
@@ -227,7 +237,6 @@ def getStartOrder(self):
         """
         return [], []
 
-
     def start(self):
         """
         This method instructs the device to begin execution of the task.
@@ -346,7 +355,7 @@ def listSequence(self):
         """
         return {}
 
-    def generateTask(self, params: "Dict | None" = None) -> dict:
+    def generateTask(self, params: Optional[dict] = None) -> dict:
         """
         This method should convert params' index-values back into task-values, along with any default work non-sequenced
         tasks need. WARNING! Long sequences will not automatically lock the UI or preserve the state of your parameter

acq4/devices/MicroManagerCamera/mmcamera.py

@@ -94,12 +94,10 @@ def _acquireFrames(self, n=1):
         self.mmc.setCameraDevice(self.camName)
         self.mmc.startSequenceAcquisition(n, 0, True)
         frames = []
-        frameTimes = []
         timeoutStart = ptime.time()
         while self.mmc.isSequenceRunning() or self.mmc.getRemainingImageCount() > 0:
             if self.mmc.getRemainingImageCount() > 0:
-                frameTimes.append(ptime.time())
-                timeoutStart = frameTimes[-1]
+                timeoutStart = ptime.time()
                 frames.append(self.mmc.popNextImage().T[np.newaxis, ...])
             else:
                 if ptime.time() - timeoutStart > 10.0:

acq4/devices/MicroManagerStage/mmstage.py

@@ -10,6 +10,7 @@
 from acq4.util.Thread import Thread
 from acq4.util.micromanager import getMMCorePy
 from ..Stage import Stage, MoveFuture, StageInterface
+from ...util.debug import printExc
 
 
 class MicroManagerStage(Stage):
@@ -131,6 +132,10 @@ def abort(self):
             except:
                 printExc("Error stopping axis %s:" % ax)
 
+    @property
+    def positionUpdatesPerSecond(self):
+        return 1 / self.monitor.minInterval
+
     def setUserSpeed(self, v):
         """Set the maximum speed of the stage (m/sec) when under manual control.
 
@@ -205,6 +210,7 @@ def __init__(self, dev):
         self.lock = Mutex(recursive=True)
         self.stopped = False
         self.interval = 0.3
+        self.minInterval = 100e-3
 
         Thread.__init__(self)
 
@@ -221,8 +227,7 @@ def setInterval(self, i):
             self.interval = i
 
     def run(self):
-        minInterval = 100e-3
-        interval = minInterval
+        interval = self.minInterval
         lastPos = None
         while True:
             try:
@@ -234,7 +239,7 @@ def run(self):
                 pos = self.dev._getPosition()  # this causes sigPositionChanged to be emitted
                 if pos != lastPos:
                     # if there was a change, then loop more rapidly for a short time.
-                    interval = minInterval
+                    interval = self.minInterval
                     lastPos = pos
                 else:
                     interval = min(maxInterval, interval * 2)

acq4/devices/Microscope/Microscope.py

@@ -1,11 +1,12 @@
-# -*- coding: utf-8 -*-
-from __future__ import print_function
-
 import collections
+from typing import Tuple
 
 import numpy as np
+import scipy.ndimage
 
 import pyqtgraph as pg
+from acq4.util.typing import Number
+from pyqtgraph.units import µm
 
 from acq4.Manager import getManager
 from acq4.devices.Device import Device
@@ -15,7 +16,7 @@
 from acq4.util import Qt
 from acq4.util.Mutex import Mutex
 from acq4.util.debug import printExc
-from acq4.util.future import MultiFuture
+from acq4.util.future import Future, MultiFuture
 
 Ui_Form = Qt.importTemplate('.deviceTemplate')
 
@@ -211,7 +212,64 @@ def setFocusDepth(self, z, speed='fast'):
         fdpos[2] += dif
         return fd.moveToGlobal(fdpos, speed)
 
-    def getSurfaceDepth(self):
+    def getDefaultImager(self):
+        name = self.config.get('defaultImager', None)
+        if name is None:
+            cameras = self.dm.listInterfaces("camera")
+            if len(cameras) == 0:
+                raise RuntimeError("No camera devices available.")
+            name = cameras[0]
+        return self.dm.getDevice(name)
+
+    def getZStack(self, imager: "Device", z_range: tuple[Number]) -> Future:
+        """Acquire a z-stack of images using the given imager.
+
+        The z-stack is returned as frames.
+        """
+        from acq4.util.imaging.sequencer import runZStack
+
+        return runZStack(imager, z_range)
+
+    @Future.wrap
+    def findSurfaceDepth(self, imager: "Device", _future: Future) -> None:
+        """Set the surface of the sample based on how focused the images are."""
+
+        def center_area(img: np.ndarray) -> Tuple[slice, slice]:
+            """Return a slice that selects the center of the image."""
+            minimum = 50
+            center_w = img.shape[0] // 2
+            start_w = max(min(int(img.shape[0] * 0.4), center_w - minimum), 0)
+            end_w = max(min(int(img.shape[0] * 0.6), center_w + minimum), img.shape[0])
+            center_h = img.shape[1] // 2
+            start_h = max(min(int(img.shape[1] * 0.4), center_h - minimum), 0)
+            end_h = max(min(int(img.shape[1] * 0.6), center_h + minimum), img.shape[1])
+            return slice(start_w, end_w), slice(start_h, end_h)
+
+        def downsample(arr, n):
+            new_shape = n * (np.array(arr.shape[1:]) / n).astype(int)
+            clipped = arr[:, :new_shape[0], :new_shape[1]]
+            mean1 = clipped.reshape(clipped.shape[0], clipped.shape[1], clipped.shape[2]//n, n).mean(axis=3)
+            mean2 = mean1.reshape(mean1.shape[0], mean1.shape[1]//n, n, mean1.shape[2]).mean(axis=2)
+            return mean2
+
+        def calculate_focus_score(image):
+            # image += np.random.normal(size=image.shape, scale=100)
+            image = scipy.ndimage.laplace(image) / np.mean(image)
+            return image.var()
+
+        z_range = (self.getSurfaceDepth() + 200 * µm, max(0, self.getSurfaceDepth() - 200 * µm), 1 * µm)
+        z_stack = _future.waitFor(self.getZStack(imager, z_range)).getResult()
+        filtered = downsample(np.array([f.data() for f in z_stack]), 5)
+        centers = filtered[(..., *center_area(filtered[0]))]
+        scored = np.array([calculate_focus_score(img) for img in centers])
+        surface = np.argmax(scored > 0.005)  # arbitrary threshold? seems about right on the test data
+        if surface == 0:
+            return
+
+        surface_frame = z_stack[surface]
+        self.setSurfaceDepth(surface_frame.info()['Depth'])
+
+    def getSurfaceDepth(self) -> Number:
         """Return the z-position of the sample surface as marked by the user.
         """
         return self._surfaceDepth
@@ -340,7 +398,6 @@ def __repr__(self):
         return "<Objective %s.%s offset=%0.2g,%0.2g scale=%0.2g>" % (self._scope.name(), self.name(), self.offset().x(), self.offset().y(), self.scale().x())
 
 
-
 class ScopeGUI(Qt.QWidget):
     """Microscope GUI displayed in Manager window.
     Shows selection of objectives and allows scale/offset to be changed for each."""
@@ -487,8 +544,12 @@ def __init__(self, dev, mod):
         self.layout.addWidget(self.setSurfaceBtn, 0, 0)
         self.setSurfaceBtn.clicked.connect(self.setSurfaceClicked)
 
+        self.findSurfaceBtn = Qt.QPushButton('Find Surface')
+        self.layout.addWidget(self.findSurfaceBtn, 1, 0)
+        self.findSurfaceBtn.clicked.connect(self.findSurfaceClicked)
+
         self.depthLabel = pg.ValueLabel(suffix='m', siPrefix=True)
-        self.layout.addWidget(self.depthLabel, 1, 0)
+        self.layout.addWidget(self.depthLabel, 2, 0)
 
         dev.sigGlobalTransformChanged.connect(self.transformChanged)
         dev.sigSurfaceDepthChanged.connect(self.surfaceDepthChanged)
@@ -504,6 +565,9 @@ def setSurfaceClicked(self):
         self.getDevice().setSurfaceDepth(focus)
         self.transformChanged()
 
+    def findSurfaceClicked(self):
+        self.getDevice().findSurfaceDepth(self.getDevice().getDefaultImager())
+
     def surfaceDepthChanged(self, depth):
         self.surfaceLine.setValue(depth)
 

acq4/devices/MockCamera/mock_camera.py

@@ -1,11 +1,9 @@
 from collections import OrderedDict
+from typing import Union
 
 import numpy as np
 import scipy
-import six
 import time
-from six.moves import range
-from six.moves import zip
 
 import acq4.util.functions as fn
 import acq4.util.ptime as ptime
@@ -127,7 +125,7 @@ def startCamera(self):
         self.lastFrameTime = ptime.time()
 
     def stopCamera(self):
-        pass
+        self.lastFrameTime = None
 
     def getNoise(self, shape):
         n = shape[0] * shape[1]
@@ -221,10 +219,30 @@ def pixelVectors(self):
 
         return x, y
 
+    def _acquireFrames(self, n: int):
+        self.startCamera()
+        try:
+            frames = []
+            while True:
+                frames.extend([f["data"][np.newaxis, ...] for f in self.newFrames()])
+                if len(frames) >= n or not self._cameraRunning():
+                    break
+                time.sleep(0.1)
+        finally:
+            self.stopCamera()
+
+        return np.concatenate(frames[:int(n)])
+
+    def _cameraRunning(self):
+        return self.lastFrameTime is not None
+
     def newFrames(self):
         """Return a list of all frames acquired since the last call to newFrames."""
         prof = pg.debug.Profiler(disabled=True)
 
+        if self.lastFrameTime is None:
+            return []
+
         now = ptime.time()
         dt = now - self.lastFrameTime
         exp = self.getParam("exposure")
@@ -284,27 +302,21 @@ def newFrames(self):
         prof()
 
         self.frameId += 1
-        frames = []
-        for i in range(nf):
-            frames.append({"data": data, "time": now + (i / fps), "id": self.frameId})
+        frames = [{"data": data, "time": now + (i / fps), "id": self.frameId} for i in range(nf)]
         prof()
         return frames
 
     def quit(self):
         pass
 
-    def listParams(self, params=None):
-        """List properties of specified parameters, or of all parameters if None"""
+    def listParams(self, params: Union[list, str, None] = None):
+        """List properties of specified parameter(s), or of all parameters if None"""
         if params is None:
             return self.paramRanges
-        else:
-            if isinstance(params, six.string_types):
-                return self.paramRanges[params]
+        if isinstance(params, str):
+            return self.paramRanges[params]
 
-            out = OrderedDict()
-            for k in params:
-                out[k] = self.paramRanges[k]
-            return out
+        return {k: self.paramRanges[k] for k in params}
 
     def setParams(self, params, autoRestart=True, autoCorrect=True):
         dp = []
@@ -363,14 +375,23 @@ def makeExpWave(self):
         cmd = self.parentTask().tasks[daq].cmd
         start = self.parentTask().startTime
         sampleRate = cmd["rate"]
-
-        data = np.zeros(cmd["numPts"], dtype=np.uint8)
-        for f in self.frames:
-            t = f.info()["time"]
-            exp = f.info()["exposure"]
-            i0 = int((t - start) * sampleRate)
-            i1 = i0 + int((exp - 0.1e-3) * sampleRate)
-            data[i0:i1] = 1
+        numPts = cmd["numPts"]
+        data = np.zeros(numPts, dtype=np.uint8)
+        if self.fixedFrameCount is None:
+            frames = self._future.peekAtResult()  # not exact, but close enough for a mock
+            for f in frames:
+                t = f.info()["time"]
+                exp = f.info()["exposure"]
+                i0 = int((t - start) * sampleRate)
+                i1 = i0 + int((exp - 0.1e-3) * sampleRate)
+                data[i0:i1] = 1
+        else:
+            n = self.fixedFrameCount
+            exp = int((self.dev.getParam("exposure") - 0.1e-3) * sampleRate)
+            minLength = max(numPts - exp, exp * n)
+            for i0 in np.linspace(1, minLength - 2, n, dtype=int):
+                i1 = i0 + exp
+                data[i0:i1] = 1
 
         return data