diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml index 4ff4e93..32eef13 100644 --- a/.github/workflows/pytest.yml +++ b/.github/workflows/pytest.yml @@ -1,7 +1,7 @@ name: PyTest on: -# [workflow_call, push, pull_request] + [workflow_call, push, pull_request] jobs: build: diff --git a/STYLEGUIDE.md b/STYLEGUIDE.md index 0f814e0..62ef254 100644 --- a/STYLEGUIDE.md +++ b/STYLEGUIDE.md @@ -6,8 +6,18 @@ # General -- 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. +- 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 -- Tests must *not* plot figures. \ No newline at end of file +- Tests must *not* plot figures. + +# Sphinx + +Common warnings: + +- All line numbers are relative to the start of the docstring. +- 'WARNING: Block quote ends without a blank line; unexpected unindent'. This happens if a block of text is not properly + wrapped and one of the lines starts with a space. To fox, remove the space at the beginning of the line. +- 'ERROR: Unexpected indentation' can be caused if a line ends with ':' and the next line is not indented or empty. \ No newline at end of file diff --git a/docs/source/references.bib b/docs/source/references.bib index dcf3513..404a9f5 100644 --- a/docs/source/references.bib +++ b/docs/source/references.bib @@ -525,11 +525,6 @@ @article{Park2012 } -@misc{pydevice, - title = {{PyDevice} {GitHub} repository}, - url = {https://www.github.com/IvoVellekoop/pydevice}, -} - @article{Pinkard2021, author = {Henry Pinkard and Nico Stuurman and Ivan E Ivanov and Nicholas M Anthony and Wei Ouyang and Bin Li and Bin Yang and Mark A Tsuchida and Bryant Chhun and Grace Zhang and Ryan Mei and Michael Anderson and Douglas P Shepherd and Ian Hunt-Isaak and Raymond L Dunn and Wiebke Jahr and Saul Kato and Loïc A Royer and Jay R Thiagarajah and Kevin W Eliceiri and Emma Lundberg and Shalin B Mehta and Laura Waller}, journal = {Nature Methods}, diff --git a/examples/micro_manager_microscope.py b/examples/micro_manager_microscope.py index 0c80dfc..2dae31c 100644 --- a/examples/micro_manager_microscope.py +++ b/examples/micro_manager_microscope.py @@ -1,5 +1,5 @@ """ Micro-Manager simulated microscope -======================= +====================================================================== This script simulates a microscope with a random noise image as a mock specimen. The numerical aperture, stage position, and other parameters can be modified through the Micro-Manager GUI. To use this script as a device in Micro-Manager, make sure you have the PyDevice adapter installed and diff --git a/examples/micro_manager_scanning_microscope.py b/examples/micro_manager_scanning_microscope.py index 5e69694..637804e 100644 --- a/examples/micro_manager_scanning_microscope.py +++ b/examples/micro_manager_scanning_microscope.py @@ -1,5 +1,5 @@ """ Micro-Manager simulated scanning microscope -======================= +====================================================================== This script simulates a scanning microscope with a pre-set image as a mock specimen. The scan parameters can be modified through the Micro-Manager GUI. To use this script as a device in Micro-Manager, make sure you have the PyDevice adapter installed and @@ -7,11 +7,11 @@ """ import astropy.units as u +import skimage # add 'openwfs' to the search path. This is only needed when developing openwfs # otherwise it is just installed as a package import set_path # noqa -import skimage from openwfs.devices import ScanningMicroscope, Axis from openwfs.devices.galvo_scanner import InputChannel diff --git a/openwfs/algorithms/basic_fourier.py b/openwfs/algorithms/basic_fourier.py index ec1d8da..822bc0e 100644 --- a/openwfs/algorithms/basic_fourier.py +++ b/openwfs/algorithms/basic_fourier.py @@ -9,14 +9,14 @@ class FourierDualReference(DualReference): - """ - Fourier double reference algorithm, based on Mastiani et al. [1]. + """Fourier double reference algorithm, based on Mastiani et al. [1]. Improvements over [1]: - * The set of plane waves is taken from a disk in k-space instead of a square. - * No overlap between the two halves is needed, instead the final stitching step is done using measurements already in the data set. - * When only a single target is optimized, the algorithm can be used in an iterative version to increase SNR during the measurument, - similar to [2]. + + - The set of plane waves is taken from a disk in k-space instead of a square. + - No overlap between the two halves is needed, instead the final stitching step is done using measurements already in the data set. + - When only a single target is optimized, the algorithm can be used in an iterative version to increase SNR during the measurument, + similar to [2]. [1]: Bahareh Mastiani, Gerwin Osnabrugge, and Ivo M. Vellekoop, "Wavefront shaping for forward scattering," Opt. Express 30, 37436-37445 (2022) diff --git a/openwfs/algorithms/custom_iter_dual_reference.py b/openwfs/algorithms/custom_iter_dual_reference.py deleted file mode 100644 index 8539289..0000000 --- a/openwfs/algorithms/custom_iter_dual_reference.py +++ /dev/null @@ -1,250 +0,0 @@ -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, - 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 b81f828..bffcff9 100644 --- a/openwfs/algorithms/dual_reference.py +++ b/openwfs/algorithms/dual_reference.py @@ -1,4 +1,4 @@ -from typing import Optional, Union +from typing import Optional, Union, List import numpy as np from numpy import ndarray as nd @@ -8,28 +8,32 @@ class DualReference: - """ - A generic iterative dual reference WFS algorithm, which can use a custom set of basis functions. + """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. + This algorithm is adapted from [Tao2017]_, 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. + 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]. + pattern as a reference. This makes this algorithm suitable for non-linear feedback, such as multi-photon + excitation fluorescence [Osnabrugge2019]_. 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). + References + ---------- + .. [Tao2017] 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). + + .. [Osnabrugge2019] 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 - [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__( @@ -148,7 +152,7 @@ def phase_patterns(self, value): b0_index = next(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: - raise ("For multi-target optimization, the both sets must contain a flat wavefront with phase 0.") + raise "For multi-target optimization, the both sets must contain a flat wavefront with phase 0." 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.") @@ -184,7 +188,7 @@ def _compute_cobasis(self): 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.") + raise "The phase_patterns must be set before computing the cobasis." cobasis = [None, None] for side in range(2): @@ -218,7 +222,7 @@ def execute(self, *, capture_intermediate_results: bool = False, progress_bar=No ref_phases = np.zeros(self._shape) # Initialize storage lists - results_all = [None] * self.iterations # List to store all results + results_all: List[Optional[WFSResult]] = [None] * self.iterations # List to store all results intermediate_results = np.zeros(self.iterations) # List to store feedback from full patterns # Prepare progress bar diff --git a/openwfs/algorithms/genetic.py b/openwfs/algorithms/genetic.py index da91588..c59b677 100644 --- a/openwfs/algorithms/genetic.py +++ b/openwfs/algorithms/genetic.py @@ -13,21 +13,28 @@ class SimpleGenetic: in [1] and [2]. The algorithm performs the following steps: + 1. Initialize all wavefronts in the population with random phases - 2. For each generation - 1. Determine the feedback signal for each wavefront - 2. Select the 'elite_size' best wavefronts to keep. Replace the rest: - * If the elite wavefronts are too similar (> 97% identical elements), - randomly generate the new wavefronts - * Otherwise, generate new wavefronts by randomly selecting two elite wavefronts - and mixing them randomly (element-wise). - Then perform a mutation that replaces a fraction (`mutation_probability`) of the - elements by a new value. - - [1]: Conkey D B, Brown A N, Caravaca-Aguirre A M and Piestun R 'Genetic algorithm optimization - for focusing through turbid media in noisy environments' Opt. Express 20 4840–9 (2012). - [2]: Benjamin R Anderson et al 'A modular GUI-based program for genetic algorithm-based - feedback-assisted wavefront shaping', J. Phys. Photonics 6 045008 (2024). + + 2. For each generation: + + 2.1. Determine the feedback signal for each wavefront. + + 2.2. Select the 'elite_size' best wavefronts to keep. Replace the rest with new wavefronts. + 2.2.1 If the elite wavefronts are too similar (> 97% identical elements), + randomly generate the new wavefronts. + + 2.2.2 Otherwise, generate new wavefronts by randomly selecting two elite wavefronts + and mixing them randomly (element-wise). + Then perform a mutation that replaces a fraction (`mutation_probability`) of the + elements by a new value. + + References + ---------- + [^1]: Conkey D B, Brown A N, Caravaca-Aguirre A M and Piestun R 'Genetic algorithm optimization + for focusing through turbid media in noisy environments' Opt. Express 20 4840–9 (2012). + [^2]: Benjamin R Anderson et al 'A modular GUI-based program for genetic algorithm-based + feedback-assisted wavefront shaping', J. Phys. Photonics 6 045008 (2024). """ def __init__( @@ -51,6 +58,7 @@ def __init__( generations (int): The number of generations mutation_probability (int): Fraction of elements in the offspring to mutate generator: a `np.random.Generator`, defaults to np.random.default_rng() + """ if np.prod(feedback.data_shape) != 1: raise ValueError("Only scalar feedback is supported") diff --git a/openwfs/algorithms/ssa.py b/openwfs/algorithms/ssa.py index 48368af..c2cebf0 100644 --- a/openwfs/algorithms/ssa.py +++ b/openwfs/algorithms/ssa.py @@ -11,7 +11,7 @@ class StepwiseSequential: TODO: enable low-res pre-optimization (as in Vellekoop & Mosk 2007) TODO: modulate segments in pupil (circle) only - [^1]:Vellekoop, I. M., & Mosk, A. P. (2007). Focusing coherent light through opaque strongly scattering media. + [1]: Vellekoop, I. M., & Mosk, A. P. (2007). Focusing coherent light through opaque strongly scattering media. Optics Letters, 32(16), 2309-2311. [2]: Ivo M. Vellekoop, "Feedback-based wavefront shaping," Opt. Express 23, 12189-12206 (2015) """ diff --git a/openwfs/devices/galvo_scanner.py b/openwfs/devices/galvo_scanner.py index 21b2259..1fe513f 100644 --- a/openwfs/devices/galvo_scanner.py +++ b/openwfs/devices/galvo_scanner.py @@ -92,18 +92,20 @@ def maximum_scan_speed(self, linear_range: float): It is assumed that the mirror accelerates and decelerates at the maximum acceleration, and scans with a constant velocity over the linear range. There are two limits to the scan speed: - * A practical limit: if it takes longer to perform the acceleration + deceleration than - it does to traverse the linear range, it does not make sense to set the scan speed so high. - The speed at which acceleration + deceleration takes as long as the linear range is the maximum speed. - * A hardware limit: when accelerating with the maximum acceleration over a distance - 0.5 * (V_max-V_min) * (1-linear_range), - the mirror will reach the maximum possible speed. + + - A practical limit: if it takes longer to perform the acceleration + deceleration than + it does to traverse the linear range, it does not make sense to set the scan speed so high. + The speed at which acceleration + deceleration takes as long as the linear range is the maximum speed. + - A hardware limit: when accelerating with the maximum acceleration over a distance + 0.5 · (V_max-V_min) · (1-linear_range), + the mirror will reach the maximum possible speed. Args: linear_range (float): fraction of the full range that is used for the linear part of the scan Returns: Quantity[u.V / u.s]: maximum scan speed + """ # x = 0.5 · a · t² = 0.5 (v_max - v_min) · (1 - linear_range) t_accel = np.sqrt((self.v_max - self.v_min) * (1 - linear_range) / self.maximum_acceleration)