Merge pull request #33 from IvoVellekoop/slm-transform

fix: SLM transpose did not take into account translation

docs/source/conf.py

@@ -55,7 +55,7 @@
         \affil[1]{University of Twente, Biomedical Photonic Imaging, TechMed Institute, P. O. Box 217,
          7500 AE Enschede, The Netherlands}
         \affil[2]{Currently at: The Netherlands Cancer Institute, Division of Molecular Pathology, 1066 CX Amsterdam, The Netherlands}
-        \affil[3]{Imec (Netherlands), Holst Centre (HTC-31), 5656 AE, Eindhoven, The Netherlands}
+        \affil[3]{Currently at: Imec (Netherlands), Holst Centre (HTC-31), 5656 AE, Eindhoven, The Netherlands}
         \affil[*]{Corresponding author: [email protected]}
         \publishers{%
             \normalfont\normalsize%
@@ -122,7 +122,7 @@
 napoleon_use_rtype = False
 napoleon_use_param = True
 typehints_document_rtype = False
-latex_engine = "pdflatex"
+latex_engine = "xelatex"
 html_theme = "sphinx_rtd_theme"
 add_module_names = False
 autodoc_preserve_defaults = True
@@ -155,7 +155,7 @@ def setup(app):
     # monkey-patch the MarkdownTranslator class to support citations
     # TODO: this should be done in the markdown builder itself
     cls = MarkdownBuilder.default_translator_class
-    cls.visit_citation = cls.visit_footnote
+    # cls.visit_citation = cls.visit_footnote
 
 
 def source_read(app, docname, source):

docs/source/readme.rst

@@ -59,7 +59,7 @@ To use OpenWFS, Python 3.9 or later is required. Since it is available on the Py
 
     pip install openwfs[all]
 
-This will also install the optional dependencies for OpenWFS:
+It is advised to make a new environment for OpenWFS, such as with conda, poetry, or Python's venv. This will also install the optional dependencies for OpenWFS:
 
 *opengl* For the OpenGL-accelerated SLM control, the ``PyOpenGL`` package is installed. In order for this package to work, an OpenGL-compatible graphics card and driver is required.
 

examples/wfs_demonstration_experimental.py

@@ -2,8 +2,14 @@
 WFS demo experiment
 =====================
 This script demonstrates how to perform a wavefront shaping experiment using the openwfs library.
-It assumes that you have a genicam camera and an SLM connected to your computer.
-Please adjust the path to the camera driver and (when needed) the monitor id in the Camera and SLM objects.
+It performs wavefront shaping with a FourierDualReference algorithm, maximizing feedback from a
+single region of interest (ROI) on a camera.
+It assumes that you have a genicam-compatible camera and an SLM connected to the video output of your computer.
+
+Adjustments to make for your specific setup:
+* the path to the camera driver.
+* The code assumes the SLM is connected as a secondary monitor. If not, adjust the monitor_id below.
+* The gray value that corresponds to a 2pi phase shift on the SLM. This is used to set the lookup table of the SLM.
 """
 
 import astropy.units as u
@@ -13,19 +19,30 @@
 from openwfs.algorithms import FourierDualReference
 from openwfs.devices import Camera, SLM
 from openwfs.processors import SingleRoi
+from openwfs.utilities import Transform
+
+# Adjust these parameters to your setup
+# The camera driver file path
+camera_driver_path = R"C:\Program Files\Basler\pylon 7\Runtime\x64\ProducerU3V.cti"
+monitor_id = 2
+# time to stabilize the SLM, in frames
+settle_time = 2
+# we are using a setup with an SLM that produces 2pi phase shift
+# at a gray value of 142. With proper calibration, this is usually 256.
+two_pi_gray_value = 142
+centering = Transform(source_origin=(0, 0), destination_origin=(0.0, 0.15))
+
 
 # This script shows how a wavefront shaping experiment can be performed from Python
-cam = Camera(R"C:\Program Files\Basler\pylon 7\Runtime\x64\ProducerU3V.cti")
+cam = Camera(camera_driver_path)
 cam.exposure_time = 16.666 * u.ms
 roi_detector = SingleRoi(cam, radius=2)
 
 # constructs the actual slm for wavefront shaping, and a monitor window to display the current phase pattern
-slm = SLM(monitor_id=2, duration=2)
-monitor = slm.clone(monitor_id=0, pos=(0, 0), shape=(slm.shape[0] // 4, slm.shape[1] // 4))
+slm = SLM(monitor_id=monitor_id, duration=settle_time, transform=centering)
+monitor = slm.clone(monitor_id=0, pos=(0, 0), shape=(slm.shape[0] // 3, slm.shape[1] // 3))
 
-# we are using a setup with an SLM that produces 2pi phase shift
-# at a gray value of 142
-slm.lookup_table = range(142)
+slm.lookup_table = range(two_pi_gray_value)
 alg = FourierDualReference(feedback=roi_detector, slm=slm, slm_shape=[800, 800], k_radius=7)
 
 result = alg.execute()

openwfs/devices/slm/slm.py

@@ -1,3 +1,5 @@
+import os
+import signal
 import warnings
 from typing import Union, Optional, Sequence
 from weakref import WeakSet
@@ -39,6 +41,7 @@ class SLM(Actuator, PhaseSLM):
         "_window",
         "_globals",
         "_frame_buffer",
+        "_monitor",
         "patches",
         "primary_patch",
         "_coordinate_system",
@@ -103,17 +106,18 @@ def __init__(
         # construct window for displaying the SLM pattern
         SLM._init_glfw()
         self._assert_window_available(monitor_id)
-        self._monitor_id = monitor_id
         self._position = pos
-        (default_shape, default_rate, _) = SLM._current_mode(monitor_id)
+        self._monitor_id = monitor_id
+        (default_shape, default_rate, _) = SLM._current_mode(self._monitor_id)
         self._shape = default_shape if shape is None else shape
         self._refresh_rate = default_rate if refresh_rate is None else refresh_rate.to_value(u.Hz)
         self._frame_buffer = None
+        self._monitor = None
         self._window = None
         self._globals = -1
         self.patches = []
         self._context = None
-        self._create_window()  # sets self._context, self._window and self._globals and self._frame_patch
+        self._create_window()  # sets self._context, self._window and self._globals and self._frame_patch, self._monitor
         self._coordinate_system = coordinate_system
         self.transform = Transform() if transform is None else transform
         self._vertex_array = VertexArray()
@@ -255,6 +259,31 @@ def _init_glfw():
         glfw.window_hint(glfw.BLUE_BITS, 8)
         glfw.window_hint(glfw.COCOA_RETINA_FRAMEBUFFER, glfw.FALSE)  # disable retina multisampling on Mac (untested)
         glfw.window_hint(glfw.SAMPLES, 0)  # disable multisampling
+        glfw.set_monitor_callback(SLM._monitor_callback)
+
+    @staticmethod
+    def _monitor_callback(monitor, event):
+        """Callback function that is called when a monitor is connected or disconnected.
+
+        If a monitor is disconnencted while an SLM window is being shown on it, the window is closed, causing all
+        subsequent access to it to fail.
+        """
+        print("Monitor event", monitor, event)
+        if event == glfw.DISCONNECTED:
+            monitor_id = glfw.get_monitor_user_pointer(monitor)
+            for slm in SLM._active_slms:
+                if slm.monitor_id == monitor_id:
+                    glfw.set_window_should_close(slm._window, 1)
+
+    @staticmethod
+    def _key_callback(window, key, scancode, action, mods):
+        """Callback function that is called when a key is pressed or released.
+
+        This callback responds to Ctrl-C by terminating the program.
+        """
+        if key == glfw.KEY_C and action == glfw.PRESS and mods == glfw.MOD_CONTROL:
+            print("Ctrl + C pressed, terminating program")
+            os.kill(os.getpid(), signal.SIGINT)
 
     def _create_window(self):
         """Constructs a new window and associated OpenGL context. Called by SLM.__init__()"""
@@ -265,12 +294,18 @@ def _create_window(self):
         shared = other._window if other is not None else None  # noqa: ok to use _window
         SLM._active_slms.add(self)
 
-        monitor = glfw.get_monitors()[self._monitor_id - 1] if self._monitor_id != SLM.WINDOWED else None
+        if self._monitor_id == SLM.WINDOWED:
+            self._monitor = None
+        else:
+            self._monitor = glfw.get_monitors()[self._monitor_id - 1]
+            glfw.set_monitor_user_pointer(self._monitor, self._monitor_id)
+
         glfw.window_hint(glfw.REFRESH_RATE, int(self._refresh_rate))
-        self._window = glfw.create_window(self._shape[1], self._shape[0], "OpenWFS SLM", monitor, shared)
+        self._window = glfw.create_window(self._shape[1], self._shape[0], "OpenWFS SLM", self._monitor, shared)
+        glfw.set_key_callback(self._window, SLM._key_callback)
         glfw.set_input_mode(self._window, glfw.CURSOR, glfw.CURSOR_HIDDEN)  # disable cursor
-        if monitor:  # full screen mode
-            glfw.set_gamma(monitor, 1.0)
+        if self._monitor:  # full screen mode
+            glfw.set_gamma(self._monitor, 1.0)
         else:  # windowed mode
             glfw.set_window_pos(self._window, self._position[1], self._position[0])
 
@@ -360,12 +395,12 @@ def monitor_id(self, value):
         self._assert_window_available(value)
         self._monitor_id = value
         (self._shape, self._refresh_rate, _) = SLM._current_mode(value)
-        monitor = glfw.get_monitors()[value - 1] if value != SLM.WINDOWED else None
+        self._monitor = glfw.get_monitors()[value - 1] if value != SLM.WINDOWED else None
 
         # move window to new monitor
         glfw.set_window_monitor(
             self._window,
-            monitor,
+            self._monitor,
             self._position[1],
             self._position[0],
             self._shape[1],
@@ -376,7 +411,8 @@ def monitor_id(self, value):
 
     def __del__(self):
         """Destructor for the SLM object. This function destroys the window and releases all resources."""
-        glfw.destroy_window(self._window)
+        if hasattr(self, "_window") and self._window is not None:  # may still be missing if window creation failed
+            glfw.destroy_window(self._window)
 
     def update(self):
         """Sends the new phase pattern to be displayed on the SLM.
@@ -411,6 +447,10 @@ def update(self):
             self._frame_buffer._draw()  # noqa - ok to access 'friend class'
 
             glfw.poll_events()  # process window messages
+            if glfw.window_should_close(self._window):
+                glfw.destroy_window(self._window)
+                self._window = None
+                raise RuntimeError("SLM window was closed")
 
             if len(self._clones) > 0:
                 self._context.__exit__(None, None, None)  # release context before updating clones

openwfs/utilities/utilities.py

@@ -42,7 +42,7 @@ def unitless(data: ArrayLike) -> np.ndarray:
 
     Usage:
     >>> data = np.array([1.0, 2.0, 3.0]) * u.m
-    >>> unitless_data = unitless(data)
+    >>> unitless_data = unitless(data / u.mm)
     """
     if isinstance(data, Quantity):
         return data.to_value(u.dimensionless_unscaled)
@@ -71,14 +71,14 @@ class Transform:
             This matrix may include astropy units,
             for instance when converting from normalized SLM coordinates to physical coordinates in micrometers.
             Note:
-                the units of the transform must match the units of the pixel_size of the destination
+                When specified, the units of the transform must match the units of the pixel_size of the destination
                 divided by the pixel size of the source.
         source_origin:
             (y,x) coordinate of the origin in the source image, relative to the center of the image.
             By default, the center of the source image is taken as the origin of the transform,
             meaning that that point is mapped onto the destination_origin.
             Note:
-                the units of source_origin must match the units of the pixel_size of the source.
+                When specified, the units of source_origin must match the units of the pixel_size of the source.
         destination_origin:
             (y,x) coordinate of the origin in the destination image, relative to the center of the image.
             By default, the center of the destination image is taken as the origin of the transform,
@@ -95,13 +95,12 @@ def __init__(
         source_origin: Optional[CoordinateType] = None,
         destination_origin: Optional[CoordinateType] = None,
     ):
-
         self.transform = Quantity(transform if transform is not None else np.eye(2))
         self.source_origin = Quantity(source_origin) if source_origin is not None else None
         self.destination_origin = Quantity(destination_origin) if destination_origin is not None else None
 
         if source_origin is not None:
-            self.destination_unit(self.source_origin.unit)  # check if the units are consistent
+            self.destination_unit(self.source_origin.unit)  # check if the units are compatible
 
     def destination_unit(self, src_unit: u.Unit) -> u.Unit:
         """Computes the unit of the output of the transformation, given the unit of the input.
@@ -126,11 +125,15 @@ def cv2_matrix(
         destination_shape: Sequence[int],
         destination_pixel_size: CoordinateType,
     ) -> np.ndarray:
-        """Returns the transformation matrix in the format used by cv2.warpAffine."""
+        """Returns the transformation matrix in the format used by cv2.warpAffine.
+
+        This matrix is 2x3, with the last column corresponding to the translation vector.
 
-        # correct the origin. OpenCV uses the _center_ of the top-left corner as the origin
-        # and openwfs uses the center of the image as the origin, so we need to shift the origin
-        # by half the image size minus half a pixel.
+        Note: OpenCV uses the (x,y) convention, so the matrix is transposed compared to the numpy matrix.
+        Note: OpenCV uses the _center_ of the top-left corner as the origin, whereas OpenWFS uses the centerso the origin is corrected.
+            of the image as the origin. The difference of half the image size minus half the
+            pixel size is corrected in the transformation matrix.
+        """
         if min(source_shape) < 1 or min(destination_shape) < 1:
             raise ValueError("Image size must be positive")
 
@@ -155,12 +158,13 @@ def cv2_matrix(
         )
 
         # finally, convert the matrix to the format used by cv2.warpAffine by swapping x and y columns and rows
-        transform_matrix = transform_matrix[[1, 0], :]
-        transform_matrix = transform_matrix[:, [1, 0, 2]]
+        transform_matrix = transform_matrix[[1, 0, 2], 0:3]
+        transform_matrix = transform_matrix[0:2, [1, 0, 2]]  # discard last row
         return transform_matrix
 
     def to_matrix(self, source_pixel_size: CoordinateType, destination_pixel_size: CoordinateType) -> np.ndarray:
-        matrix = np.zeros((2, 3))
+        """Returns a homogeneous transformation matrix that transforms (y,x,1) coordinates to (y', x', 1)."""
+        matrix = np.eye(3)
         matrix[0:2, 0:2] = unitless(self.transform * source_pixel_size / destination_pixel_size)
         if self.destination_origin is not None:
             matrix[0:2, 2] = unitless(self.destination_origin / destination_pixel_size)
@@ -169,22 +173,22 @@ def to_matrix(self, source_pixel_size: CoordinateType, destination_pixel_size: C
         return matrix
 
     def opencl_matrix(self) -> np.ndarray:
-        # compute the transform for points (0,1), (1,0), and (0, 0)
-        matrix = self.to_matrix((1.0, 1.0), (1.0, 1.0))
+        """Returns the transformation matrix (including translation) in the format used by OpenCL.
 
-        # subtract the offset (transform of 0,0) from the first two points
-        # to construct the homogeneous transformation matrix
-        # convert to opencl format: swap x and y columns (note: the rows were
-        # already swapped in the construction of t2), and flip the sign of the y-axis.
+        The returned matrix is a 3x4 matrix that includes the transformation matrix and the translation vector.
+        The last two columns are for padding only.
+        """
+
+        # Compute the homogeneous transform matrix, and transpose it.
+        # Then, swap x and y rows and columns,
+        # and flip the sign of the y-axis
+        matrix = self.to_matrix((1.0, 1.0), (1.0, 1.0)).transpose()
+
+        # convert to opencl format: swap x and y rows and columns,
+        # and flip the sign of the y-axis output.
         transform = np.eye(3, 4, dtype="float32", order="C")
-        transform[0, 0:3] = matrix[
-            1,
-            [1, 0, 2],
-        ]
-        transform[1, 0:3] = -matrix[
-            0,
-            [1, 0, 2],
-        ]
+        transform[:, 0] = matrix[[1, 0, 2], 1]
+        transform[:, 1] = -matrix[[1, 0, 2], 0]
         return transform
 
     @staticmethod
@@ -240,9 +244,10 @@ def compose(self, other: "Transform"):
             Transform: the composition of the two transformations
         """
         transform = self.transform @ other.transform
-        source_origin = other.source_origin
-        destination_origin = self.apply(other.destination_origin) if other.destination_origin is not None else None
-        return Transform(transform, source_origin, destination_origin)
+        destination_origin = (
+            self.apply(other.destination_origin) if other.destination_origin is not None else self.destination_origin
+        )
+        return Transform(transform, other.source_origin, destination_origin)
 
     def _standard_input(self) -> Quantity:
         """Construct standard input points (1,0), (0,1) and (0,0) with the source unit of this transform."""

tests/test_utilities.py

@@ -19,13 +19,13 @@ def test_to_matrix():
     )
 
     # Define the expected output matrix for same input and output pixel sizes
-    expected_matrix = ((1, 2, 1), (3, 4, 1))
+    expected_matrix = ((1, 2, 1), (3, 4, 1), (0, 0, 1))
     result_matrix = transform.to_matrix((1, 2) * u.um, (1, 2) * u.um)
-    assert result_matrix.shape == (2, 3)
+    assert result_matrix.shape == (3, 3)
     assert np.allclose(result_matrix, expected_matrix)
 
     # Repeat for different input and output pixel sizes
-    expected_matrix = ((0.5, 4, 1), (1.5, 8, 1))
+    expected_matrix = ((0.5, 4, 1), (1.5, 8, 1), (0, 0, 1))
     result_matrix = transform.to_matrix((0.5, 4) * u.um, (1, 2) * u.um)
     assert np.allclose(result_matrix, expected_matrix)
 
@@ -54,13 +54,13 @@ def test_to_matrix():
     assert np.allclose(result_matrix @ src_center, dst_center)
 
     # Also check openGL matrix (has y-axis flipped and extra row and column)
-    expected_matrix = ((1, 2, 1), (3, 4, 2))
+    expected_matrix = ((1, 2, 1), (3, 4, 2), (0, 0, 1))
     transform = Transform(transform=((1, 2), (3, 4)), source_origin=(0, 0), destination_origin=(1, 2))
     result_matrix = transform.to_matrix((1, 1), (1, 1))
     assert np.allclose(result_matrix, expected_matrix)
 
     result_matrix = transform.opencl_matrix()
-    expected_matrix = ((4, 3, 2, 0), (-2, -1, -1, 0), (0, 0, 1, 0))
+    expected_matrix = ((4, -2, 0, 0), (3, -1, 0, 0), (2, -1, 1, 0))
     assert np.allclose(result_matrix, expected_matrix)