diff --git a/.github/workflows/black.yml b/.github/workflows/black.yml new file mode 100644 index 0000000..68b72e8 --- /dev/null +++ b/.github/workflows/black.yml @@ -0,0 +1,10 @@ +name: Black + +on: [workflow_call, push, pull_request] + +jobs: + lint: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: psf/black@stable \ No newline at end of file diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml new file mode 100644 index 0000000..d0b9728 --- /dev/null +++ b/.github/workflows/pytest.yml @@ -0,0 +1,24 @@ +name: PyTest + +on: [workflow_call, push, pull_request] + +jobs: + build: + runs-on: ubuntu-latest + strategy: + matrix: + python-version: ["3.9", "3.10", "3.11"] + + steps: + - uses: actions/checkout@v4 + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install .[dev] + - name: Test with pytest + run: | + pytest \ No newline at end of file diff --git a/README.md b/README.md index 2051751..8409041 100644 --- a/README.md +++ b/README.md @@ -4,6 +4,11 @@ +[![PyTest]([https://github.com/IvoVellekoop/openwfs/actions/workflows/pytest.yml/badge.svg)](https://github.com/IvoVellekoop/openwfs/actions/workflows/pytest.yml](https://github.com/IvoVellekoop/openwfs/actions/workflows/pytest.yml/badge.svg)](https://github.com/IvoVellekoop/openwfs/actions/workflows/pytest.yml)) +[![Black]([https://github.com/IvoVellekoop/openwfs/actions/workflows/black.yml/badge.svg)](https://github.com/IvoVellekoop/openwfs/actions/workflows/black.yml](https://github.com/IvoVellekoop/openwfs/actions/workflows/black.yml/badge.svg)](https://github.com/IvoVellekoop/openwfs/actions/workflows/black.yml)) + +What is wavefront shaping? + ## What is wavefront shaping? Wavefront shaping (WFS) is a technique for controlling the propagation of light in arbitrarily complex structures, including strongly scattering materials [[1](#id62)]. In WFS, a spatial light modulator (SLM) is used to shape the phase and/or amplitude of the incident light. With a properly constructed wavefront, light can be made to focus through [[2](#id48)], or inside [[3](#id37)] scattering materials; or light can be shaped to have other desired properties, such as optimal sensitivity for specific measurements [[4](#id38)], specialized point-spread functions [[5](#id25)] or for functions like optical trapping [[6](#id28)]. diff --git a/STYLEGUIDE.md b/STYLEGUIDE.md index ed33cef..0f814e0 100644 --- a/STYLEGUIDE.md +++ b/STYLEGUIDE.md @@ -6,9 +6,7 @@ # General -- all .py files MUST be formatted with the 'black' autoformatter. This can be done by installing the 'black' package, - and running `black .` in the root directory. black is automatically installed when the development dependencies are - included. +- The package `black` is used to ensure correct formatting. Install with `pip install black` and run in the terminal using `black .` when located at the root of the repository. # Tests diff --git a/docs/source/development.rst b/docs/source/development.rst index f2b1a1c..7b57fcb 100644 --- a/docs/source/development.rst +++ b/docs/source/development.rst @@ -14,7 +14,9 @@ To download the source code, including tests and examples, clone the repository poetry install --with dev --with docs poetry run pytest -The examples are located in the ``examples`` directory. Note that a lot of functionality is also demonstrated in the automatic tests located in the ``tests`` directory. As an alternative to downloading the source code, the samples can also be copied directly from the example gallery on the documentation website :cite:`readthedocsOpenWFS`. +The examples are located in the ``examples`` directory. Note that a lot of functionality is also demonstrated in the automatic tests located in the ``tests`` directory. As an alternative to downloading the source code, the samples can also be copied directly from the example gallery on the documentation website :cite:`readthedocsOpenWFS`. + +Important to note for adding hardware devices, is that many of the components rely on third-party, in some case proprietary drivers. For using NI DAQ components, the nidaqmx package needs to be installed, and for openCV and Genicam their respective drivers need to be installed. The specific requirements are always listed in the documentation of the functions and classes that require packages like these. Building the documentation -------------------------------------------------- diff --git a/docs/source/readme.rst b/docs/source/readme.rst index ebd0fad..65cf83d 100644 --- a/docs/source/readme.rst +++ b/docs/source/readme.rst @@ -12,15 +12,19 @@ OpenWFS :target: https://openwfs.readthedocs.io/en/latest/?badge=latest :alt: Documentation Status +.. only:: markdown + + [![PyTest](https://github.com/IvoVellekoop/openwfs/actions/workflows/pytest.yml/badge.svg)](https://github.com/IvoVellekoop/openwfs/actions/workflows/pytest.yml) + [![Black](https://github.com/IvoVellekoop/openwfs/actions/workflows/black.yml/badge.svg)](https://github.com/IvoVellekoop/openwfs/actions/workflows/black.yml) What is wavefront shaping? -------------------------------- -Wavefront shaping (WFS) is a technique for controlling the propagation of light in arbitrarily complex structures, including strongly scattering materials :cite:`kubby2019`. In WFS, a spatial light modulator (SLM) is used to shape the phase and/or amplitude of the incident light. With a properly constructed wavefront, light can be made to focus through :cite:`Vellekoop2007`, or inside :cite:`vellekoop2008demixing` scattering materials; or light can be shaped to have other desired properties, such as optimal sensitivity for specific measurements :cite:`bouchet2021maximum`, specialized point-spread functions :cite:`boniface2017transmission` or for functions like optical trapping :cite:`vcivzmar2010situ`. +Wavefront shaping (WFS) is a technique for controlling the propagation of light in arbitrarily complex structures, including strongly scattering materials :cite:`kubby2019`. In WFS, a spatial light modulator (SLM) is used to shape the phase and/or amplitude of the incident light. With a properly constructed wavefront, light can be made to focus through :cite:`Vellekoop2007`, or inside :cite:`vellekoop2008demixing` scattering materials; or light can be shaped to have other desired properties, such as optimal sensitivity for specific measurements :cite:`bouchet2021maximum`, specialized point-spread functions :cite:`boniface2017transmission`, spectral filtering :cite:`Park2012`,, or for functions like optical trapping :cite:`vcivzmar2010situ`. -It stands out that an important driving force in WFS is the development of new algorithms, for example to account for sample movement :cite:`valzania2023online`, to be optimally resilient to noise :cite:`mastiani2021noise`, or to use digital twin models to compute the required correction patterns :cite:`salter2014exploring,ploschner2015seeing,Thendiyammal2020,cox2023model`. Much progress has been made towards developing fast and noise-resilient algorithms, or algorithms designed for specific towards the methodology of wavefront shaping, such as using algorithms based on Hadamard patterns, or Fourier-based approaches :cite:`Mastiani2022`. Fast techniques that enable wavefront shaping in dynamic samples :cite:`Liu2017,Tzang2019`, and many potential applications have been developed and prototyped, including endoscopy :cite:`ploschner2015seeing`, optical trapping :cite:`Cizmar2010` and deep-tissue imaging :cite:`Streich2021`. +It stands out that an important driving force in WFS is the development of new algorithms, for example to account for sample movement :cite:`valzania2023online`, experimental conditions :cite:`Anderson2016`, to be optimally resilient to noise :cite:`mastiani2021noise`, or to use digital twin models to compute the required correction patterns :cite:`salter2014exploring,ploschner2015seeing,Thendiyammal2020,cox2023model`. Much progress has been made towards developing fast and noise-resilient algorithms, or algorithms designed for specific towards the methodology of wavefront shaping, such as using algorithms based on Hadamard patterns, or Fourier-based approaches :cite:`Mastiani2022`. Fast techniques that enable wavefront shaping in dynamic samples :cite:`Liu2017,Tzang2019`, and many potential applications have been developed and prototyped, including endoscopy :cite:`ploschner2015seeing`, optical trapping :cite:`Cizmar2010`, Raman scattering, :cite:`Thompson2016`, and deep-tissue imaging :cite:`Streich2021`. Applications extend beyond that of microscope imaging such as optimizing photoelectrochemical absorption :cite:`Liew2016` and tuning random lasers :cite:`Bachelard2014`. -With the development of these advanced algorithms, however, the complexity of WFS software is gradually becoming a bottleneck for further advancements in the field, as well as for end-user adoption. Code for controlling wavefront shaping tends to be complex and setup-specific, and developing this code typically requires detailed technical knowledge and low-level programming. Moreover, since many labs use their own in-house programs to control the experiments, sharing and re-using code between different research groups is troublesome. +With the development of these advanced algorithms, however, the complexity of WFS software is steadily increasing as the field matures, which hinders cooperation as well as end-user adoption. Code for controlling wavefront shaping tends to be complex and setup-specific, and developing this code typically requires detailed technical knowledge and low-level programming. A recent c++ based contribution :cite:`Anderson2024`, highlights the growing need for software based tools that enable use and development. Moreover, since many labs use their own in-house programs to control the experiments, sharing and re-using code between different research groups is troublesome. What is OpenWFS? ---------------------- @@ -34,11 +38,13 @@ OpenWFS is a Python package for performing and for simulating wavefront shaping * **GenICam cameras**. The :class:`~.devices.Camera` object uses the `harvesters` backend :cite:`harvesters` to access any camera supporting the GenICam standard :cite:`genicam`. * **Automatic synchronization**. OpenWFS provides tools for automatic synchronization of actuators (e. g. an SLM) and detectors (e. g. a camera). The automatic synchronization makes it trivial to perform pipelined measurements that avoid the delay normally caused by the latency of the video card and SLM. -* **Wavefront shaping algorithms**. A (growing) collection of wavefront shaping algorithms. OpenWFS abstracts the hardware control, synchronization, and signal processing so that the user can focus on the algorithm itself. As a result, most algorithms can be implemented in just a few lines of code without the need for low-level or hardware-specific programming. +* **Wavefront shaping algorithms**. A (growing) collection of wavefront shaping algorithms. OpenWFS abstracts the hardware control, synchronization, and signal processing so that the user can focus on the algorithm itself. As a result, most algorithms can be implemented cleanly without hardware-specific programming. * **Simulation**. OpenWFS provides an extensive framework for testing and simulating wavefront shaping algorithms, including the effect of measurement noise, stage drift, and user-defined aberrations. This allows for rapid prototyping and testing of new algorithms, without the need for physical hardware. -* **Platform for exchange and joint collaboration**. OpenWFS can be used as a platform for sharing and exchanging wavefront shaping algorithms. The package is designed to be modular and easy to expand, and it is our hope that the community will contribute to the package by adding new algorithms, hardware control modules, and simulation tools. +* **Platform for exchange and joint collaboration**. OpenWFS can be used as a platform for sharing and exchanging wavefront shaping algorithms. The package is designed to be modular and easy to expand, and it is our hope that the community will contribute to the package by adding new algorithms, hardware control modules, and simulation tools. Python was specifically chosen for this purpose for its active community, high level of abstraction and the ease of sharing tools. Further expansion of the supported hardware is of high priority, especially wrapping c-based software support with tools like ctypes and the Micro-Manager based device adapters. + +* **Platform for simplifying use of wavefront shaping**. OpenWFS is compatible with the recently developed PyDevice :cite:`PyDevice`, and can therefore be controlled from Micro-Manager :cite:`MMoverview`, a commonly used microscopy control platform. * **Automated troubleshooting**. OpenWFS provides tools for automated troubleshooting of wavefront shaping experiments. This includes tools for measuring the performance of wavefront shaping algorithms, and for identifying common problems such as incorrect SLM calibration, drift, measurement noise, and other experimental imperfections. diff --git a/docs/source/references.bib b/docs/source/references.bib index fe4be20..998292c 100644 --- a/docs/source/references.bib +++ b/docs/source/references.bib @@ -157,6 +157,24 @@ @article{astropy primaryClass = {astro-ph.IM}, } +@article{Thompson2016, + abstract = {Spontaneous Raman scattering is a powerful tool for chemical sensing and imaging but suffers from a weak signal. In this Letter, we present an application of adaptive optics to enhance the Raman scattering signal detected through a turbid, optically thick material. This technique utilizes recent advances in wavefront shaping techniques for focusing light through a turbid media and applies them to chemical detection to achieve a signal enhancement with little sacrifice to the overall simplicity of the experimental setup. With this technique, we demonstrate an enhancement in the Raman signal from titanium dioxide particles through a highly scattering material. This technique may pave the way to label-free tracking using the optical memory effect.}, + author = {Jonathan V. Thompson and Graham A. Throckmorton and Brett H. Hokr and Vladislav V. Yakovlev}, + doi = {10.1364/OL.41.001769}, + issn = {1539-4794}, + issue = {8}, + journal = {Optics letters}, + keywords = {Graham A Throckmorton,Jonathan V Thompson,MEDLINE,NCBI,NIH,NLM,National Center for Biotechnology Information,National Institutes of Health,National Library of Medicine,Non-P.H.S.,PubMed Abstract,Radiation*,Raman*,Research Support,Scattering,Spectrum Analysis,U.S. Gov't,Vladislav V Yakovlev,doi:10.1364/OL.41.001769,pmid:27082341}, + month = {4}, + pages = {1769}, + pmid = {27082341}, + publisher = {Opt Lett}, + title = {Wavefront shaping enhanced Raman scattering in a turbid medium}, + volume = {41}, + url = {https://pubmed.ncbi.nlm.nih.gov/27082341/}, + year = {2016}, +} + @article{vellekoop2008demixing, title = {Demixing light paths inside disordered metamaterials}, @@ -364,6 +382,41 @@ @article{Streich2021 year = {2021}, } +@article{Liew2016, + abstract = {A fundamental issue that limits the efficiency of many photoelectrochemical systems is that the photon absorption length is typically much longer than the electron diffusion length. Various photon management schemes have been developed to enhance light absorption; one simple approach is to use randomly scattering media to enable broadband and wide-angle enhancement. However, such systems are often opaque, making it difficult to probe photoinduced processes. Here we use wave interference effects to modify the spatial distribution of light inside a highly scattering dye-sensitized solar cell to control photon absorption in a space-dependent manner. By shaping the incident wavefront of a laser beam, we enhance or suppress photocurrent by increasing or decreasing light concentration on the front side of the mesoporous photoanode where the collection efficiency of photoelectrons is maximal. Enhanced light absorption is achieved by reducing reflection through the open boundary of the photoanode via destructive interference, leading to a factor of 2 increase in photocurrent. This approach opens the door to probing and manipulating photoelectrochemical processes in specific regions inside nominally opaque media.}, + author = {Seng Fatt Liew and Sébastien M. Popoff and Stafford W. Sheehan and Arthur Goetschy and Charles A. Schmuttenmaer and A. Douglas Stone and Hui Cao}, + doi = {10.1021/ACSPHOTONICS.5B00642}, + issn = {2330-4022}, + issue = {3}, + journal = {ACS Photonics}, + keywords = {dye-sensitized solar cells,multimode interference,multiple scattering,photoelectrochemical,wavefront shaping}, + month = {3}, + pages = {449-455}, + publisher = {American Chemical Society}, + title = {Coherent Control of Photocurrent in a Strongly Scattering Photoelectrochemical System}, + volume = {3}, + url = {https://technion-staging.elsevierpure.com/en/publications/coherent-control-of-photocurrent-in-a-strongly-scattering-photoel}, + year = {2016}, +} + + +@article{Anderson2016, + abstract = {Previously we considered the effect of experimental parameters on optimized transmission through opaque media using spatial light modulator (SLM)-based wavefront shaping. In this study we consider the opposite geometry, in which we optimize reflection from an opaque surface such that the backscattered light is focused onto a spot on an imaging detector. By systematically varying different experimental parameters (genetic algorithm iterations, bin size, SLM active area, target area, spot size, and sample angle with respect to the optical axis) and optimizing the reflected light we determine how each parameter affects the intensity enhancement. We find that the effects of the experimental parameters on the enhancement are similar to those measured for a transmissive geometry, but with the exact functional forms changed due to the different geometry and the use of a genetic algorithm instead of an iterative algorithm. Additionally, we find preliminary evidence of greater enhancements than predicted by random matrix theory, suggesting a possibly new physical mechanism to be investigated in future work.}, + author = {Benjamin R. Anderson and Ray Gunawidjaja and Hergen Eilers}, + doi = {10.1103/PHYSREVA.93.013813/FIGURES/12/MEDIUM}, + issn = {24699934}, + issue = {1}, + journal = {Physical Review A}, + month = {1}, + pages = {013813}, + publisher = {American Physical Society}, + title = {Effect of experimental parameters on optimal reflection of light from opaque media}, + volume = {93}, + url = {https://journals.aps.org/pra/abstract/10.1103/PhysRevA.93.013813}, + year = {2016}, +} + + @misc{MMoverview, author = {Mark Tsuchida and Sam Griffin}, title = {Micro-Manager Project Overview}, @@ -375,6 +428,60 @@ @misc{openwfsdocumentation url = {https://openwfs.readthedocs.io/en/latest/}, } +@article{Anderson2024, + abstract = {We have developed a modular graphical user interface (GUI)-based program for use in genetic algorithm-based feedback-assisted wavefront shaping. The program uses a class-based structure to separate out the universal modules (e.g. GUI, multithreading, optimization algorithms) and hardware-specific modules (e.g. code for different SLMs and cameras). This modular design makes the program easily adaptable to a wide range of lab equipment, while providing easy access to a GUI, multithreading, and three optimization algorithms (phase-stepping, simple genetic, and microgenetic).}, + author = {Benjamin R. Anderson and Andrew O’Kins and Kostiantyn Makrasnov and Rebecca Udby and Patrick Price and Hergen Eilers}, + doi = {10.1088/2515-7647/AD6ED3}, + issn = {2515-7647}, + issue = {4}, + journal = {Journal of Physics: Photonics}, + keywords = {disordered media,genetic algorithms,multithreading,spatial light modulators,wavefront shaping}, + month = {8}, + pages = {045008}, + publisher = {IOP Publishing}, + title = {A modular GUI-based program for genetic algorithm-based feedback-assisted wavefront shaping}, + volume = {6}, + url = {https://iopscience.iop.org/article/10.1088/2515-7647/ad6ed3 https://iopscience.iop.org/article/10.1088/2515-7647/ad6ed3/meta}, + year = {2024}, +} + + +@article{Bachelard2014, + abstract = {A laser is not necessarily a sophisticated device: pumping an amplifying medium randomly filled with scatterers makes a perfectly viable â ̃ random laserâ ™. The absence of mirrors greatly simplifies laser design, but control over the emission wavelength and directionality is lost, seriously hindering prospects for this otherwise simple laser. Recently, we proposed an approach to tame random lasers, inspired by coherent light control in complex media. Here, we implement this method in an optofluidic random laser where modes are spatially extended and overlap, making individual mode selection impossible, a priori. We show experimentally that control over laser emission can be regained even in this extreme case. By actively shaping the optical pump within the random laser, single-mode operation at any selected wavelength is achieved with spectral selectivity down to 0.06 nm and more than 10 dB side-lobe rejection. This method paves the way towards versatile tunable and controlled random lasers as well as the taming of other laser sources. © 2014 Macmillan Publishers Limited. All rights reserved.}, + author = {Nicolas Bachelard and Sylvain Gigan and Xavier Noblin and Patrick Sebbah}, + doi = {10.1038/nphys2939}, + issn = {17452481}, + issue = {6}, + journal = {Nature Physics}, + keywords = {Optics,Physics}, + pages = {426-431}, + publisher = {Nature Publishing Group}, + title = {Adaptive pumping for spectral control of random lasers}, + volume = {10}, + url = {https://ui.adsabs.harvard.edu/abs/2014NatPh..10..426B/abstract}, + year = {2014}, +} + + +@article{Park2012, + abstract = {We demonstrate controlled wavelength-dependent light focusing through turbid media using wavefront shaping. Due to the dispersion caused by multiple light scattering, light propagation through turbid media can be independently controlled between different wavelengths. Foci with various wavelengths can be generated by applying different optimized wavefronts to a highly scattering layer. Given the linearity of the transmission matrix, multiple foci with different wavelengths can also be simultaneously constructed.}, + author = {Jung-Hoon Park and ChungHyun Park and YongKeun Park and Hyunseung Yu and Yong-Hoon Cho}, + doi = {10.1364/OL.37.003261}, + issn = {1539-4794}, + issue = {15}, + journal = {Optics Letters, Vol. 37, Issue 15, pp. 3261-3263}, + keywords = {Light propagation,Multiple scattering,Scattering media,Second harmonic generation,Spatial light modulators,Turbid media}, + month = {8}, + pages = {3261-3263}, + pmid = {22859152}, + publisher = {Optica Publishing Group}, + title = {Active spectral filtering through turbid media}, + volume = {37}, + url = {https://opg.optica.org/viewmedia.cfm?uri=ol-37-15-3261&seq=0&html=true https://opg.optica.org/abstract.cfm?uri=ol-37-15-3261 https://opg.optica.org/ol/abstract.cfm?uri=ol-37-15-3261}, + year = {2012}, +} + + @misc{pydevice, title = {{PyDevice} {GitHub} repository}, url = {https://www.github.com/IvoVellekoop/pydevice}, @@ -452,3 +559,10 @@ @book{kubby2019 editor = {Kubby, Joel and Gigan, Sylvain and Cui, Meng}, year = {2019}, collection = {Advances in Microscopy and Microanalysis} } + +@misc{PyDevice, + author = {Ivo Vellekoop and Jeroen Doornbos}, + title = {Micro-Manager PyDevice}, + url = {https://micro-manager.org/PyDevice}, +} + diff --git a/openwfs/algorithms/custom_iter_dual_reference.py b/openwfs/algorithms/custom_iter_dual_reference.py new file mode 100644 index 0000000..73b083b --- /dev/null +++ b/openwfs/algorithms/custom_iter_dual_reference.py @@ -0,0 +1,220 @@ +from typing import Optional + +import numpy as np +from numpy import ndarray as nd + +from .utilities import analyze_phase_stepping, WFSResult +from ..core import Detector, PhaseSLM + + +def weighted_average(a, b, wa, wb): + """ + Compute the weighted average of two values. + + Args: + a: The first value. + b: The second value. + wa: The weight of the first value. + wb: The weight of the second value. + """ + return (a * wa + b * wb) / (wa + wb) + + +class IterativeDualReference: + """ + A generic iterative dual reference WFS algorithm, which can use a custom set of basis functions. + + This algorithm is adapted from [1], with the addition of the ability to use custom basis functions and specify the number of iterations. + + In this algorithm, the SLM pixels are divided into two groups: A and B, as indicated by the boolean group_mask argument. + The algorithm first keeps the pixels in group B fixed, and displays a sequence on patterns on the pixels of group A. + It uses these measurements to construct an optimized wavefront that is displayed on the pixels of group A. + This process is then repeated for the pixels of group B, now using the *optimized* wavefront on group A as reference. + Optionally, the process can be repeated for a number of iterations, which each iteration using the current correction + pattern as a reference. This makes this algorithm suitable for non-linear feedback, such as multi-photon + excitation fluorescence [2]. + + This algorithm assumes a phase-only SLM. Hence, the input modes are defined by passing the corresponding phase + patterns (in radians) as input argument. + + [1]: X. Tao, T. Lam, B. Zhu, et al., “Three-dimensional focusing through scattering media using conjugate adaptive + optics with remote focusing (CAORF),” Opt. Express 25, 10368–10383 (2017). + + [2]: Gerwin Osnabrugge, Lyubov V. Amitonova, and Ivo M. Vellekoop. "Blind focusing through strongly scattering media + using wavefront shaping with nonlinear feedback", Optics Express, 27(8):11673–11688, 2019. + https://opg.optica.org/oe/ abstract.cfm?uri=oe-27-8-1167 + """ + + def __init__(self, feedback: Detector, slm: PhaseSLM, phase_patterns: tuple[nd, nd], group_mask: nd, + phase_steps: int = 4, iterations: int = 4, analyzer: Optional[callable] = analyze_phase_stepping): + """ + Args: + feedback: The feedback source, usually a detector that provides measurement data. + slm: Spatial light modulator object. + phase_patterns: A tuple of two 3D arrays, containing the phase patterns for group A and group B, respectively. + The first two dimensions are the spatial dimensions, and should match the size of group_mask. + The 3rd dimension in the array is index of the phase pattern. The number of phase patterns in A and B may be different. + group_mask: A 2D bool array of that defines the pixels used by group A with False and elements used by + group B with True. + phase_steps: The number of phase steps for each mode (default is 4). Depending on the type of + non-linear feedback and the SNR, more might be required. + iterations: Number of times to measure a mode set, e.g. when iterations = 5, the measurements are + A, B, A, B, A. Should be at least 2 + analyzer: The function used to analyze the phase stepping data. Must return a WFSResult object. Defaults to `analyze_phase_stepping` + """ + if (phase_patterns[0].shape[0:2] != group_mask.shape) or (phase_patterns[1].shape[0:2] != group_mask.shape): + raise ValueError("The phase patterns and group mask must all have the same shape.") + if iterations < 2: + raise ValueError("The number of iterations must be at least 2.") + if np.prod(feedback.data_shape) != 1: + raise ValueError("The feedback detector should return a single scalar value.") + + self.slm = slm + self.feedback = feedback + self.phase_steps = phase_steps + self.iterations = iterations + self.analyzer = analyzer + self.phase_patterns = (phase_patterns[0].astype(np.float32), phase_patterns[1].astype(np.float32)) + mask = group_mask.astype(bool) + self.masks = (~mask, mask) # masks[0] is True for group A, mask[1] is True for group B + + # Pre-compute the conjugate modes for reconstruction + self.modes = [np.exp(-1j * self.phase_patterns[side]) * np.expand_dims(self.masks[side], axis=2) for side in + range(2)] + + def execute(self, capture_intermediate_results: bool = False, progress_bar=None) -> WFSResult: + """ + Executes the blind focusing dual reference algorithm and compute the SLM transmission matrix. + capture_intermediate_results: When True, measures the feedback from the optimized wavefront after each iteration. + This can be useful to determine how many iterations are needed to converge to an optimal pattern. + This data is stored as the 'intermediate_results' field in the results + progress_bar: Optional progress bar object. Following the convention for tqdm progress bars, + this object should have a `total` attribute and an `update()` function. + + Returns: + WFSResult: An object containing the computed SLM transmission matrix and related data. The amplitude profile + of each mode is assumed to be 1. If a different amplitude profile is desired, this can be obtained by + multiplying that amplitude profile with this transmission matrix. + """ + + # Current estimate of the transmission matrix (start with all 0) + t_full = np.zeros(shape=self.modes[0].shape[0:2]) + t_other_side = t_full + + # Initialize storage lists + t_set_all = [None] * self.iterations + results_all = [None] * self.iterations # List to store all results + results_latest = [None, None] # The two latest results. Used for computing fidelity factors. + intermediate_results = np.zeros(self.iterations) # List to store feedback from full patterns + + # Prepare progress bar + if progress_bar: + num_measurements = np.ceil(self.iterations / 2) * self.modes[0].shape[2] \ + + np.floor(self.iterations / 2) * self.modes[1].shape[2] + progress_bar.total = num_measurements + + # Switch the phase sets back and forth multiple times + for it in range(self.iterations): + side = it % 2 # pick set A or B for phase stepping + ref_phases = -np.angle(t_full) # use the best estimate so far to construct an optimized reference + side_mask = self.masks[side] + # Perform WFS experiment on one side, keeping the other side sized at the ref_phases + result = self._single_side_experiment(mod_phases=self.phase_patterns[side], ref_phases=ref_phases, + mod_mask=side_mask, progress_bar=progress_bar) + + # Compute transmission matrix for the current side and update + # estimated transmission matrix + t_this_side = self.compute_t_set(result, self.modes[side]) + t_full = t_this_side + t_other_side + t_other_side = t_this_side + + # Store results + t_set_all[it] = t_this_side # Store transmission matrix + results_all[it] = result # Store result + results_latest[side] = result # Store latest result for this set + + # Try full pattern + if capture_intermediate_results: + self.slm.set_phases(-np.angle(t_full)) + intermediate_results[it] = self.feedback.read() + + # Compute average fidelity factors + fidelity_noise = weighted_average(results_latest[0].fidelity_noise, + results_latest[1].fidelity_noise, results_latest[0].n, + results_latest[1].n) + fidelity_amplitude = weighted_average(results_latest[0].fidelity_amplitude, + results_latest[1].fidelity_amplitude, results_latest[0].n, + results_latest[1].n) + fidelity_calibration = weighted_average(results_latest[0].fidelity_calibration, + results_latest[1].fidelity_calibration, results_latest[0].n, + results_latest[1].n) + + result = WFSResult(t=t_full, + t_f=None, + n=self.modes[0].shape[2] + self.modes[1].shape[2], + axis=2, + fidelity_noise=fidelity_noise, + fidelity_amplitude=fidelity_amplitude, + fidelity_calibration=fidelity_calibration) + + # TODO: document the t_set_all and results_all attributes + result.t_set_all = t_set_all + result.results_all = results_all + result.intermediate_results = intermediate_results + return result + + def _single_side_experiment(self, mod_phases: nd, ref_phases: nd, mod_mask: nd, + progress_bar=None) -> WFSResult: + """ + Conducts experiments on one part of the SLM. + + Args: + mod_phases: 3D array containing the phase patterns of each mode. Axis 0 and 1 are used as spatial axis. + Axis 2 is used for the 'phase pattern index' or 'mode index'. + ref_phases: 2D array containing the reference phase pattern. + mod_mask: 2D array containing a boolean mask, where True indicates the modulated part of the SLM. + progress_bar: Optional progress bar object. Following the convention for tqdm progress bars, + this object should have a `total` attribute and an `update()` function. + + Returns: + WFSResult: An object containing the computed SLM transmission matrix and related data. + """ + num_modes = mod_phases.shape[2] + measurements = np.zeros((num_modes, self.phase_steps)) + + for m in range(num_modes): + phases = ref_phases.copy() + modulated = mod_phases[:, :, m] + for p in range(self.phase_steps): + phi = p * 2 * np.pi / self.phase_steps + # set the modulated pixel values to the values corresponding to mode m and phase offset phi + phases[mod_mask] = modulated[mod_mask] + phi + self.slm.set_phases(phases) + self.feedback.trigger(out=measurements[m, p, ...]) + + if progress_bar is not None: + progress_bar.update() + + self.feedback.wait() + return self.analyzer(measurements, axis=1) + + @staticmethod + def compute_t_set(wfs_result: WFSResult, mode_set: nd) -> nd: + """ + Compute the transmission matrix in SLM space from transmission matrix in input mode space. + + Note 1: This function computes the transmission matrix for one mode set, and thus returns one part of the full + transmission matrix. The elements that are not part of the mode set will be 0. The full transmission matrix can + be obtained by simply adding the parts, i.e. t_full = t_set0 + t_set1. + + Note 2: As this is a blind focusing WFS algorithm, there may be only one target or 'output mode'. + + Args: + wfs_result (WFSResult): The result of the WFS algorithm. This contains the transmission matrix in the space + of input modes. + mode_set: 3D array with set of modes. + """ + t = wfs_result.t.squeeze().reshape((1, 1, mode_set.shape[2])) + norm_factor = np.prod(mode_set.shape[0:2]) + t_set = (t * mode_set).sum(axis=2) / norm_factor + return t_set diff --git a/openwfs/algorithms/dual_reference.py b/openwfs/algorithms/dual_reference.py index f2fc8f5..a16d2f7 100644 --- a/openwfs/algorithms/dual_reference.py +++ b/openwfs/algorithms/dual_reference.py @@ -38,6 +38,7 @@ def __init__( feedback: Detector, slm: PhaseSLM, phase_patterns: Optional[tuple[nd, nd]], + amplitude: Optional[tuple[nd, nd] | str], group_mask: nd, phase_steps: int = 4, iterations: int = 2, @@ -52,6 +53,11 @@ def __init__( The first two dimensions are the spatial dimensions, and should match the size of group_mask. The 3rd dimension in the array is index of the phase pattern. The number of phase patterns in A and B may be different. When None, the phase_patterns attribute must be set before executing the algorithm. + amplitude: Tuple of 2D arrays, one array for each group. The arrays have shape equal to the shape of + group_mask. When None, the amplitude attribute must be set before executing the algorithm. When + 'uniform', a 2D array of normalized uniform values is used, such that ⟨A,A⟩=1, where ⟨.,.⟩ denotes the + inner product and A is the amplitude profile per group. This corresponds to a uniform illumination of + the SLM. Note: if the groups have different sizes, their normalization factors will be different. group_mask: A 2D bool array of that defines the pixels used by group A with False and elements used by group B with True. phase_steps: The number of phase steps for each mode (default is 4). Depending on the type of @@ -72,7 +78,6 @@ def __init__( [1]: X. Tao, T. Lam, B. Zhu, et al., “Three-dimensional focusing through scattering media using conjugate adaptive optics with remote focusing (CAORF),” Opt. Express 25, 10368–10383 (2017). - """ if optimized_reference is None: # 'auto' mode optimized_reference = np.prod(feedback.data_shape) == 1 @@ -95,14 +100,40 @@ def __init__( self.iterations = iterations self._analyzer = analyzer self._phase_patterns = None + self._amplitude = None + self._gram = None self._shape = group_mask.shape mask = group_mask.astype(bool) self.masks = ( ~mask, mask, - ) # mask[0] is True for group A, mask[1] is True for group B + ) # self.masks[0] is True for group A, self.masks[1] is True for group B + self.amplitude = amplitude # Note: when 'uniform' is passed, the shape of self.masks[0] is used. self.phase_patterns = phase_patterns + @property + def amplitude(self) -> Optional[nd]: + return self._amplitude + + @amplitude.setter + def amplitude(self, value): + if value is None: + self._amplitude = None + return + + if value == 'uniform': + self._amplitude = tuple( + (np.ones(shape=self._shape) / np.sqrt(self.masks[side].sum())).astype(np.float32) for side in range(2)) + return + + if value.shape != self._shape: + raise ValueError( + "The amplitude and group mask must all have the same shape." + ) + + self._amplitude = value.astype(np.float32) + + @property def phase_patterns(self) -> tuple[nd, nd]: return self._phase_patterns @@ -118,10 +149,14 @@ def phase_patterns(self, value): # find the modes in A and B that correspond to flat wavefronts with phase 0 try: a0_index = next( - i for i in range(value[0].shape[2]) if np.allclose(value[0][:, :, i], 0) + i + for i in range(value[0].shape[2]) + if np.allclose(value[0][:, :, i], 0) ) b0_index = next( - i for i in range(value[1].shape[2]) if np.allclose(value[1][:, :, i], 0) + i + for i in range(value[1].shape[2]) + if np.allclose(value[1][:, :, i], 0) ) self.zero_indices = (a0_index, b0_index) except StopIteration: @@ -130,13 +165,56 @@ def phase_patterns(self, value): ) if (value[0].shape[0:2] != self._shape) or (value[1].shape[0:2] != self._shape): - raise ValueError("The phase patterns and group mask must all have the same shape.") + raise ValueError( + "The phase patterns and group mask must all have the same shape." + ) self._phase_patterns = ( value[0].astype(np.float32), value[1].astype(np.float32), ) + self._compute_cobasis() + + @property + def cobasis(self) -> tuple[nd, nd]: + """ + The cobasis corresponding to the given basis. + """ + return self._cobasis + + @property + def gram(self) -> np.matrix: + """ + The Gram matrix corresponding to the given basis (i.e. phase pattern and amplitude profile). + """ + return self._gram + + def _compute_cobasis(self): + """ + Computes the cobasis from the phase patterns. + + As a basis matrix is full rank, this is equivalent to the Moore-Penrose pseudo-inverse. + B⁺ = (B^* B)^(-1) B^* + Where B is the basis matrix (a column corresponds to a basis vector), ^* denotes the conjugate transpose, ^(-1) + denotes the matrix inverse, and ⁺ denotes the Moore-Penrose pseudo-inverse. + """ + if self.phase_patterns is None: + raise('The phase_patterns must be set before computing the cobasis.') + + cobasis = [None, None] + for side in range(2): + p = np.prod(self._shape) # Number of SLM pixels + m = self.phase_patterns[side].shape[2] # Number of modes + phase_factor = np.exp(1j * self.phase_patterns[side]) + amplitude_factor = np.expand_dims(self.amplitude[side] * self.masks[side], axis=2) + B = np.asmatrix((phase_factor * amplitude_factor).reshape((p, m))) # Basis matrix + self._gram = B.H @ B + B_pinv = np.linalg.inv(self.gram) @ B.H # Moore-Penrose pseudo-inverse + cobasis[side] = np.asarray(B_pinv.T).reshape(self.phase_patterns[side].shape) + + self._cobasis = cobasis + def execute( self, *, capture_intermediate_results: bool = False, progress_bar=None ) -> WFSResult: @@ -155,11 +233,6 @@ def execute( """ # Current estimate of the transmission matrix (start with all 0) - cobasis = [ - np.exp(-1j * self.phase_patterns[side]) * np.expand_dims(self.masks[side], axis=2) - for side in range(2) - ] - ref_phases = np.zeros(self._shape) # Initialize storage lists @@ -193,7 +266,7 @@ def execute( if self.optimized_reference: # use the best estimate so far to construct an optimized reference - t_this_side = self.compute_t_set(results_all[it].t, cobasis[side]).squeeze() + t_this_side = self.compute_t_set(results_all[it].t, self.cobasis[side]).squeeze() ref_phases[self.masks[side]] = -np.angle(t_this_side[self.masks[side]]) # Try full pattern @@ -211,10 +284,12 @@ def execute( relative = results_all[0].t[self.zero_indices[0], ...] + np.conjugate( results_all[1].t[self.zero_indices[1], ...] ) - factor = (relative / np.abs(relative)).reshape((1, *self.feedback.data_shape)) + factor = (relative / np.abs(relative)).reshape( + (1, *self.feedback.data_shape) + ) - t_full = self.compute_t_set(results_all[0].t, cobasis[0]) + self.compute_t_set( - factor * results_all[1].t, cobasis[1] + t_full = self.compute_t_set(results_all[0].t, self.cobasis[0]) + self.compute_t_set( + factor * results_all[1].t, self.cobasis[1] ) # Compute average fidelity factors @@ -247,7 +322,9 @@ def _single_side_experiment( WFSResult: An object containing the computed SLM transmission matrix and related data. """ num_modes = mod_phases.shape[2] - measurements = np.zeros((num_modes, self.phase_steps, *self.feedback.data_shape)) + measurements = np.zeros( + (num_modes, self.phase_steps, *self.feedback.data_shape) + ) for m in range(num_modes): phases = ref_phases.copy() diff --git a/openwfs/plot_utilities.py b/openwfs/plot_utilities.py index 1a9a6bd..e11d1e8 100644 --- a/openwfs/plot_utilities.py +++ b/openwfs/plot_utilities.py @@ -1,5 +1,11 @@ +from typing import Tuple + +import numpy as np +from numpy import ndarray as nd from astropy import units as u from matplotlib import pyplot as plt +from matplotlib.colors import hsv_to_rgb +from matplotlib.axes import Axes from .core import Detector from .utilities import get_extent @@ -48,3 +54,153 @@ def scale_prefix(value: u.Quantity) -> u.Quantity: return value.to(u.s) else: return value + + +def slope_step(a: nd, width: nd | float) -> nd: + """ + A sloped step function from 0 to 1. + + Args: + a: Input array + width: width of the sloped step. + + Returns: + An array the size of a, with the result of the sloped step function. + """ + return (a >= width) + a/width * (0 < a) * (a < width) + + +def linear_blend(a: nd, b: nd, blend: nd | float) -> nd: + """ + Return a linear, element-wise blend between two arrays a and b. + + Args: + a: Input array a. + b: Input array b. + blend: Blend factor. Value of 1.0 -> return a. Value of 0.0 -> return b. + + Returns: + A linear combination of a and b, corresponding to the blend factor. a*blend + b*(1-blend) + """ + return a*blend + b*(1-blend) + + +def complex_to_rgb(array: nd, scale: float | nd | None = None, axis: int = 2) -> nd: + """ + Generate RGB color values to represent values of a complex array. + + The complex values are mapped to HSV colorspace and then converted to RGB. Hue represents phase and Value represents + amplitude. Saturation is set to 1. + + Args: + array: Array to create RGB values for. + scale: Scaling factor for the array values. When None, scale = 1/max(abs(array)) is used. + axis: Array axis to use for the RGB dimension. + + Returns: + An RGB array representing the complex input array. + """ + if scale is None: + scale = 1 / np.max(abs(array)) + h = np.expand_dims(np.angle(array) / (2 * np.pi) + 0.5, axis=axis) + s = np.ones_like(h) + v = np.expand_dims(np.abs(array) * scale, axis=axis).clip(min=0, max=1) + hsv = np.concatenate((h, s, v), axis=axis) + rgb = hsv_to_rgb(hsv) + return rgb + + +def plot_field(array, scale: float | nd | None = None, imshow_kwargs: dict | None = None): + """ + Plot a complex array as an RGB image. + + The phase is represented by the hue, and the magnitude by the value, i.e. black = zero, brightness shows amplitude, + and the colors represent the phase. + + Args: + array(ndarray): complex array to be plotted. + scale(float): scaling factor for the magnitude. The final value is clipped to the range [0, 1]. + imshow_kwargs: Keyword arguments for matplotlib's imshow. + """ + if imshow_kwargs is None: + imshow_kwargs = {} + rgb = complex_to_rgb(array, scale) + plt.imshow(rgb, **imshow_kwargs) + + +def plot_scatter_field(x, y, array, scale, scatter_kwargs=None): + """ + Plot complex scattered data as RGB values. + """ + if scatter_kwargs is None: + scatter_kwargs = {'s': 80} + rgb = complex_to_rgb(array, scale, axis=1) + plt.scatter(x, y, c=rgb, **scatter_kwargs) + + +def complex_colorbar(scale, width_inverse: int = 15): + """ + Create an rgb colorbar for complex numbers and return its Axes handle. + """ + amp = np.linspace(0, 1.01, 10).reshape((1, -1)) + phase = np.linspace(0, 249 / 250 * 2 * np.pi, 250).reshape(-1, 1) - np.pi + z = amp * np.exp(1j * phase) + rgb = complex_to_rgb(z, 1) + ax = plt.subplot(1, width_inverse, width_inverse) + plt.imshow(rgb, aspect='auto', extent=(0, scale, -np.pi, np.pi)) + + # Ticks and labels + ax.set_yticks((-np.pi, -np.pi / 2, 0, np.pi / 2, np.pi), ('$-\\pi$', '$-\\pi/2$', '0', '$\\pi/2$', '$\\pi$')) + ax.set_xlabel('amp.') + ax.set_ylabel('phase (rad)') + ax.yaxis.tick_right() + ax.yaxis.set_label_position("right") + return ax + + +def complex_colorwheel(ax: Axes = None, shape: Tuple[int, int] = (100, 100), imshow_kwargs: dict = {}, + arrow_props: dict = {}, text_kwargs: dict = {}, amplitude_str: str = 'A', + phase_str: str = '$\\phi$'): + """ + Create an rgb image for a colorwheel representing the complex unit circle. + + Args: + ax: Matplotlib Axes. + shape: Number of pixels in each dimension. + imshow_kwargs: Keyword arguments for matplotlib's imshow. + arrow_props: Keyword arguments for the arrows. + text_kwargs: Keyword arguments for the text labels. + amplitude_str: Text label for the amplitude arrow. + phase_str: Text label for the phase arrow. + + Returns: + rgb_wheel: rgb image of the colorwheel. + """ + if ax is None: + ax = plt.gca() + + x = np.linspace(-1, 1, shape[1]).reshape(1, -1) + y = np.linspace(-1, 1, shape[0]).reshape(-1, 1) + z = x + 1j*y + rgb = complex_to_rgb(z, scale=1) + step_width = 1.5 / shape[1] + blend = np.expand_dims(slope_step(1 - np.abs(z) - step_width, width=step_width), axis=2) + rgba_wheel = np.concatenate((rgb, blend), axis=2) + ax.imshow(rgba_wheel, extent=(-1, 1, -1, 1), **imshow_kwargs) + + # Add arrows with annotations + ax.annotate('', xy=(-0.98/np.sqrt(2),)*2, xytext=(0, 0), arrowprops={'color': 'white', 'width': 1.8, + 'headwidth': 5.0, 'headlength': 6.0, **arrow_props}) + ax.text(**{'x': -0.4, 'y': -0.8, 's': amplitude_str, 'color': 'white', 'fontsize': 15, **text_kwargs}) + ax.annotate('', xy=(0, 0.9), xytext=(0.9, 0), + arrowprops={'connectionstyle': 'arc3,rad=0.4', 'color': 'white', 'width': 1.8, 'headwidth': 5.0, + 'headlength': 6.0, **arrow_props}) + ax.text(**{'x': 0.1, 'y': 0.5, 's': phase_str, 'color': 'white', 'fontsize': 15, **text_kwargs}) + + # Hide axes spines and ticks + ax.set_xticks([]) + ax.set_yticks([]) + ax.spines['left'].set_visible(False) + ax.spines['right'].set_visible(False) + ax.spines['top'].set_visible(False) + ax.spines['bottom'].set_visible(False) diff --git a/tests/test_wfs.py b/tests/test_wfs.py index d04992d..17574b9 100644 --- a/tests/test_wfs.py +++ b/tests/test_wfs.py @@ -14,8 +14,9 @@ from ..openwfs.algorithms.troubleshoot import field_correlation from ..openwfs.algorithms.utilities import WFSController from ..openwfs.processors import SingleRoi -from ..openwfs.simulation import SimulatedWFS, StaticSource, SLM, Microscope from ..openwfs.simulation.mockdevices import GaussianNoise +from ..openwfs.simulation import SimulatedWFS, StaticSource, SLM, Microscope +from ..openwfs.plot_utilities import plot_field @pytest.mark.parametrize("shape", [(4, 7), (10, 7), (20, 31)]) @@ -59,7 +60,9 @@ def test_multi_target_algorithms(shape, noise: float, algorithm: str): k_radius=(np.min(shape) - 1) // 2, phase_steps=phase_steps, ) - N = alg.phase_patterns[0].shape[2] + alg.phase_patterns[1].shape[2] # number of input modes + N = ( + alg.phase_patterns[0].shape[2] + alg.phase_patterns[1].shape[2] + ) # number of input modes alg_fidelity = 1.0 # Fourier is accurate for any N signal = 1 / 2 # for estimating SNR. @@ -81,7 +84,8 @@ def test_multi_target_algorithms(shape, noise: float, algorithm: str): I_opt[b] = feedback.read()[b] t_correlation += abs(np.vdot(result.t[:, :, b], sim.t[:, :, b])) ** 2 t_norm += abs( - np.vdot(result.t[:, :, b], result.t[:, :, b]) * np.vdot(sim.t[:, :, b], sim.t[:, :, b]) + np.vdot(result.t[:, :, b], result.t[:, :, b]) + * np.vdot(sim.t[:, :, b], sim.t[:, :, b]) ) t_correlation /= t_norm @@ -92,10 +96,14 @@ def test_multi_target_algorithms(shape, noise: float, algorithm: str): # the fidelity of the transmission matrix reconstruction theoretical_noise_fidelity = signal / (signal + noise**2 / phase_steps) enhancement = I_opt.mean() / I_0 - theoretical_enhancement = np.pi / 4 * theoretical_noise_fidelity * alg_fidelity * (N - 1) + 1 + theoretical_enhancement = ( + np.pi / 4 * theoretical_noise_fidelity * alg_fidelity * (N - 1) + 1 + ) estimated_enhancement = result.estimated_enhancement.mean() * alg_fidelity theoretical_t_correlation = theoretical_noise_fidelity * alg_fidelity - estimated_t_correlation = result.fidelity_noise * result.fidelity_calibration * alg_fidelity + estimated_t_correlation = ( + result.fidelity_noise * result.fidelity_calibration * alg_fidelity + ) tolerance = 2.0 / np.sqrt(M) print( f"\nenhancement: \ttheoretical= {theoretical_enhancement},\testimated={estimated_enhancement},\tactual: {enhancement}" @@ -164,13 +172,15 @@ def test_fourier2(): @pytest.mark.skip("Not implemented") def test_fourier_microscope(): aberration_phase = skimage.data.camera() * ((2 * np.pi) / 255.0) + np.pi - aberration = StaticSource(aberration_phase, pixel_size=2.0 / np.array(aberration_phase.shape)) + aberration = StaticSource( + aberration_phase, pixel_size=2.0 / np.array(aberration_phase.shape) + ) img = np.zeros((1000, 1000), dtype=np.int16) signal_location = (250, 250) img[signal_location] = 100 slm_shape = (1000, 1000) - src = StaticSource(img, pixel_size=400 * u.nm) + src = StaticSource(img, 400 * u.nm) slm = SLM(shape=(1000, 1000)) sim = Microscope( source=src, @@ -192,8 +202,12 @@ def test_fourier_microscope(): after = roi_detector.read() # imshow(controller._optimized_wavefront) print(after / before) - scaled_aberration = zoom(aberration_phase, np.array(slm_shape) / aberration_phase.shape) - assert_enhancement(slm, roi_detector, controller._result, np.exp(1j * scaled_aberration)) + scaled_aberration = zoom( + aberration_phase, np.array(slm_shape) / aberration_phase.shape + ) + assert_enhancement( + slm, roi_detector, controller._result, np.exp(1j * scaled_aberration) + ) def test_fourier_correction_field(): @@ -212,7 +226,9 @@ def test_fourier_correction_field(): t = alg.execute().t t_correct = np.exp(1j * aberrations) - correlation = np.vdot(t, t_correct) / np.sqrt(np.vdot(t, t) * np.vdot(t_correct, t_correct)) + correlation = np.vdot(t, t_correct) / np.sqrt( + np.vdot(t, t) * np.vdot(t_correct, t_correct) + ) # TODO: integrate with other test cases, duplication assert abs(correlation) > 0.75 @@ -385,19 +401,21 @@ def test_simple_genetic(population_size: int, elite_size: int): assert after / before > 4 -@pytest.mark.parametrize("type", ("plane_wave", "hadamard")) +@pytest.mark.parametrize("basis_str", ("plane_wave", "hadamard")) @pytest.mark.parametrize("shape", ((8, 8), (16, 4))) -def test_custom_blind_dual_reference_ortho_split(type: str, shape): +def test_custom_blind_dual_reference_ortho_split(basis_str: str, shape): """Test custom blind dual reference with an orthonormal phase-only basis. Two types of bases are tested: plane waves and Hadamard""" do_debug = False N = shape[0] * (shape[1] // 2) modes_shape = (shape[0], shape[1] // 2, N) - if type == "plane_wave": + if basis_str == "plane_wave": # Create a full plane wave basis for one half of the SLM. - modes = np.fft.fft2(np.eye(N).reshape(modes_shape), axes=(0, 1)) - else: # type == 'hadamard': - modes = hadamard(N).reshape(modes_shape) + modes = np.fft.fft2(np.eye(N).reshape(modes_shape), axes=(0, 1)) / np.sqrt(N) + elif basis_str == 'hadamard': + modes = hadamard(N).reshape(modes_shape) / np.sqrt(N) + else: + raise f'Unknown type of basis "{basis_str}".' mask = np.concatenate( (np.zeros(modes_shape[0:2], dtype=bool), np.ones(modes_shape[0:2], dtype=bool)), @@ -412,12 +430,13 @@ def test_custom_blind_dual_reference_ortho_split(type: str, shape): plt.figure(figsize=(12, 7)) for m in range(N): - plt.subplot(*modes_shape[0:1], m + 1) - plt.imshow(np.angle(mode_set[:, :, m]), vmin=-np.pi, vmax=np.pi) + plt.subplot(*modes_shape[0:2], m + 1) + plot_field(mode_set[:, :, m]) plt.title(f"m={m}") plt.xticks([]) plt.yticks([]) - plt.pause(0.1) + plt.suptitle('Basis') + plt.pause(0.01) # Create aberrations sim = SimulatedWFS(t=random_transmission_matrix(shape)) @@ -426,6 +445,7 @@ def test_custom_blind_dual_reference_ortho_split(type: str, shape): feedback=sim, slm=sim.slm, phase_patterns=(phases_set, np.flip(phases_set, axis=1)), + amplitude='uniform', group_mask=mask, iterations=4, ) @@ -433,6 +453,14 @@ def test_custom_blind_dual_reference_ortho_split(type: str, shape): result = alg.execute() if do_debug: + plt.figure() + for m in range(N): + plt.subplot(*modes_shape[0:2], m + 1) + plot_field(alg.cobasis[0][:, :, m]) + plt.title(f'{m}') + plt.suptitle('Cobasis') + plt.pause(0.01) + plt.figure() plt.imshow(np.angle(sim.t), vmin=-np.pi, vmax=np.pi, cmap="hsv") plt.title("Aberrations") @@ -443,9 +471,12 @@ def test_custom_blind_dual_reference_ortho_split(type: str, shape): plt.colorbar() plt.show() - assert ( - np.abs(field_correlation(sim.t, result.t)) > 0.99 - ) # todo: find out why this is not higher + # Checks for orthonormal bases + assert np.allclose(alg.gram, np.eye(N), atol=1e-6) # Gram matrix must be I + assert np.allclose(alg.cobasis[0], mode_set.conj(), atol=1e-6) # Cobasis vectors are just the complex conjugates + + # todo: find out why this is not higher + assert np.abs(field_correlation(sim.t, result.t)) > 0.95 def test_custom_blind_dual_reference_non_ortho(): @@ -458,7 +489,7 @@ def test_custom_blind_dual_reference_non_ortho(): N1 = 6 N2 = 3 M = N1 * N2 - mode_set_half = (1 / M) * (1j * np.eye(M).reshape((N1, N2, M)) * -np.ones(shape=(N1, N2, M))) + mode_set_half = np.exp(2j*np.pi/3 * np.eye(M).reshape((N1, N2, M))) / np.sqrt(M) mode_set = np.concatenate((mode_set_half, np.zeros(shape=(N1, N2, M))), axis=1) phases_set = np.angle(mode_set) mask = np.concatenate((np.zeros((N1, N2)), np.ones((N1, N2))), axis=1) @@ -470,8 +501,8 @@ def test_custom_blind_dual_reference_non_ortho(): plt.figure(figsize=(12, 7)) for m in range(M): plt.subplot(N2, N1, m + 1) - plt.imshow(phases_set[:, :, m], vmin=-np.pi, vmax=np.pi) - plt.title(f"m={m}") + plot_field(mode_set[:, :, m]) + plt.title(f'm={m}') plt.xticks([]) plt.yticks([]) plt.pause(0.01) @@ -481,7 +512,9 @@ def test_custom_blind_dual_reference_non_ortho(): x = np.linspace(-1, 1, 1 * N1).reshape((1, -1)) y = np.linspace(-1, 1, 1 * N1).reshape((-1, 1)) aberrations = ( - np.sin(0.8 * np.pi * x) * np.cos(1.3 * np.pi * y) * (0.8 * np.pi + 0.4 * x + 0.4 * y) + np.sin(0.8 * np.pi * x) + * np.cos(1.3 * np.pi * y) + * (0.8 * np.pi + 0.4 * x + 0.4 * y) ) % (2 * np.pi) aberrations[0:1, :] = 0 aberrations[:, 0:2] = 0 @@ -492,6 +525,7 @@ def test_custom_blind_dual_reference_non_ortho(): feedback=sim, slm=sim.slm, phase_patterns=(phases_set, np.flip(phases_set, axis=1)), + amplitude='uniform', group_mask=mask, phase_steps=4, iterations=4, @@ -499,16 +533,31 @@ def test_custom_blind_dual_reference_non_ortho(): result = alg.execute() + aberration_field = np.exp(1j * aberrations) + t_field = np.exp(1j * np.angle(result.t)) + if do_debug: plt.figure() - plt.imshow(np.angle(np.exp(1j * aberrations)), vmin=-np.pi, vmax=np.pi, cmap="hsv") - plt.title("Aberrations") - plt.colorbar() + for m in range(M): + plt.subplot(N2, N1, m + 1) + plot_field(alg.cobasis[0][:, :, m], scale=2) + plt.title(f'{m}') + plt.suptitle('Cobasis') + plt.pause(0.01) plt.figure() - plt.imshow(np.angle(result.t), vmin=-np.pi, vmax=np.pi, cmap="hsv") - plt.title("t") + plt.imshow(abs(alg.gram), vmin=0, vmax=1) + plt.title('Gram matrix abs values') plt.colorbar() + + plt.figure() + plt.subplot(1, 2, 1) + plot_field(aberration_field) + plt.title('Aberrations') + + plt.subplot(1, 2, 2) + plot_field(t_field) + plt.title('t') plt.show() - assert np.abs(field_correlation(np.exp(1j * aberrations), result.t)) > 0.999 + assert np.abs(field_correlation(aberration_field, t_field)) > 0.999