diff --git a/docs/source/api.rst b/docs/source/api.rst index 340abfc..4c5dbc0 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -65,5 +65,79 @@ The unification of these two categories is stored in: cameraslms +API Formalism +============= + +Conventions +~~~~~~~~~~~ +:mod:`slmsuite` follows the ``shape = (h, w)`` and ``vector = (x, y)`` formalism adopted by +the :mod:`numpy` ecosystem. :mod:`numpy`, :mod:`scipy`, :mod:`matplotlib`, etc generally follow this +formalism. The ``shape`` and indexing of an array or image always uses the inverted ``(h, w)`` form, +but other functions such as ``numpy.meshgrid(x, y)`` (default), ``scipy.odr.Data(x, y)``, or +``matplotlib.pyplot.scatter(x, y)`` use the standard cartesian ``(x, y)`` form that is more familiar +to users. This is not ideal and causes confusion, but this is the formalism generally +adopted by the community. + +We additionally adopt the convention that +a list of :math:`N` vectors with dimension :math:`D` is represented with shape +``(D, N)``. This is so that linear transformations can be done with a direct and easy +multiplication with a ``(D, D)`` matrix. +Stacks of :math:`W \times H` images, however, +are stored with the inverted shape ``(N, H, W)``. +These conventions are arguably consistent with the chosen ordering. + +Bases +~~~~~ +This package uses a number of spatial bases or coordinate spaces. Some coordinate spaces are +directly used by the user (most often the camera basis ``"ij"`` used for feedback). +Other bases are less often used directly, but are important to how holograms are +optimized under the hood (esp. ``"knm"``, the coordinate space of optimization when +using discrete Fourier transforms). + +.. list-table:: Bases used in :mod:`slmsuite`. + :widths: 20 80 + :header-rows: 1 + + * - Basis + - Meaning + * - ``"kxy"`` + - Normalized basis of the SLM's :math:`k`-space in normalized units. + Centered at ``(kx, ky) = (0, 0)``. This basis is what the SLM projects in angular + space (which maps to the camera's image space via the Fourier transform + implemented by free space and lens separating the two). + + The edge of Fourier space corresponds to when + :math:`|k_x|` or :math:`|k_y|` equal 1. + Note that the small angle approximation breaks down in this limit. + * - ``"knm"`` + - Nyquist basis of the SLM's computational :math:`k`-space for a given + discrete computational grid of shape ``shape``. + Centered at ``(kn, km) = (shape[1]/2, shape[0]/2)``. + ``"knm"`` is a discrete version of the continuous ``"kxy"``. This is + important because holograms need to be stored in computer memory, a discrete + medium with pixels, rather than being purely continuous. For instance, in + :class:`SpotHologram`, spots targeting specific continuous angles are rounded to + the nearest discrete pixels of ``"knm"`` space in practice + (though :class:`CompressedSpotHologram` transcends this limitation). + Then, this ``"knm"`` space image is handled as a + standard image/array, and operations such as the discrete Fourier transform + (instrumental for numerical hologram optimization) can be applied. + + The edge of ``"knm"`` space is bounded by zero and the extent of ``shape``. + In ``"kxy"`` space, this edge corresponds to the Nyquist limit of the SLM + and is strictly smaller than the full extent of Fourier space. Increasing the + ``shape`` of ``"knm"`` space increases the resolution of the grid in Fourier + space, as the edge of ``"knm"`` space is fixed by the SLM. + * - ``"ij"`` + - Pixel basis of the camera. + Centered at ``(i, j) = (cam.shape[1]/2, cam.shape[0]/2)``. + Is in the image space of the camera. + + The bounds of pixel space may be larger or smaller than Fourier or Nyquist space, + depending upon the imaging optics that separate the camera and SLM. + +See the first tip in :class:`Hologram` to learn more about ``"kxy"`` and ``"knm"`` +space. + .. |slmsuite| replace:: :mod:`slmsuite` .. _slmsuite: https://github.com/QPG-MIT/slmsuite \ No newline at end of file diff --git a/slmsuite/hardware/cameraslms.py b/slmsuite/hardware/cameraslms.py index 1f37d41..59d577d 100644 --- a/slmsuite/hardware/cameraslms.py +++ b/slmsuite/hardware/cameraslms.py @@ -84,6 +84,8 @@ def __init__(self, cam, slm, mag=1): self.mag = float(mag) + self.calibrations = {} + class NearfieldSLM(CameraSLM): """ @@ -103,33 +105,50 @@ def __init__(self, cam, slm, mag=None): class FourierSLM(CameraSLM): - """ + r""" Class for an SLM and camera separated by a Fourier transform. - This class includes methods for Fourier space and wavefront calibration. + This class includes methods for system calibration. Attributes ---------- - fourier_calibration : dict OR None - Information for the affine transformation that maps between - the k-space of the SLM (kxy) and the pixel-space of the camera (ij). - See :meth:`~slmsuite.hardware.cameraslms.FourierSLM.fourier_calibrate()`. - In the future, including pincushion or other distortions might be implemented - using :meth:`scipy.ndimage.spline_filter()`, - though this mapping would be separable. - wavefront_calibration_raw : dict OR None - Raw data for wavefront calibration, which corrects for aberrations - in the optical system (phase) and measures the amplitude distribution - of power on the SLM (amplitude). - See :meth:`~slmsuite.hardware.cameraslms.FourierSLM.wavefront_calibrate`. + calibrations : dict + "fourier" : dict + The affine transformation that maps between + the k-space of the SLM (kxy) and the pixel-space of the camera (ij). + See :meth:`~slmsuite.hardware.cameraslms.FourierSLM.fourier_calibrate()`. + + This data is critical for much of :mod:`slmsuite`'s functionality. + "wavefront" : dict + Raw data for correcting aberrations in the optical system (``phase``) and + measuring the optical amplitude distribution incident on the SLM (``amp``). + See :meth:`~slmsuite.hardware.cameraslms.FourierSLM.wavefront_calibrate()`. + Usable data is produced by running + See :meth:`~slmsuite.hardware.cameraslms.FourierSLM.process_wavefront_calibration()`. + + This data is critical for crisp holography. + "pixel" : dict + Raw data for measuring the crosstalk and :math:`V_\pi` of sections of the + SLM via a + Usable data is produced by running + See :meth:`~slmsuite.hardware.cameraslms.FourierSLM.process_pixel_calibration()`. + + **This data is currently unused; trying to figure out a + computationally-efficient way to apply the crosstalk without oversampling.** + "settle" : dict + Raw data for determining the system response of the SLM. + See :meth:`~slmsuite.hardware.cameraslms.FourierSLM.settle_calibrate()`. + Usable data is produced by running + See :meth:`~slmsuite.hardware.cameraslms.FourierSLM.process_settle_calibration()`. + + This data informs the user's choice of `settle_time_s`, or the time to wait to + acquire data after a pattern is displayed. This is, of course, a tradeoff + between measurement speed and precision. """ def __init__(self, *args, **kwargs): # """See :attr:`CameraSLM.__init__`.""" super().__init__(*args, **kwargs) - self.fourier_calibration = None - self.wavefront_calibration_raw = None - # Size of the calibration point window relative to the spot. self._wavefront_calibration_window_multiplier = 4 @@ -150,7 +169,7 @@ def simulate(self): hardware. """ # Make sure we have a Fourier calibration. - if self.fourier_calibration is None: + if self.calibrations["fourier"] is None: raise ValueError("Cannot simulate() a FourierSLM without a Fourier calibration.") # Make a simulated SLM @@ -168,8 +187,8 @@ def simulate(self): cam_sim = SimulatedCamera( slm_sim, resolution=self.cam.shape[::-1], - M=copy.copy(self.fourier_calibration["M"]), - b=copy.copy(self.fourier_calibration["b"]), + M=copy.copy(self.calibrations["fourier"]["M"]), + b=copy.copy(self.calibrations["fourier"]["b"]), bitdepth=self.cam.bitdepth, averaging=self.cam._buffer.shape[0] # TODO check this. if self.cam._buffer is not None else None, @@ -180,71 +199,199 @@ def simulate(self): #Combine the two and pass FourierSLM attributes from hardware fs_sim = FourierSLM(cam_sim, slm_sim) - fs_sim.fourier_calibration = copy.deepcopy(self.fourier_calibration) - fs_sim.wavefront_calibration_raw = copy.deepcopy(self.wavefront_calibration_raw) + fs_sim.calibrations = copy.deepcopy(self.calibrations) fs_sim._wavefront_calibration_window_multiplier = self._wavefront_calibration_window_multiplier return fs_sim + ### Calibration Helpers ### + + def name_calibration(self, calibration_type): + """ + Creates ``"--"``. + + Parameters + ---------- + calibration_type : str + The type of calibration to save. See :attr:`calibrations` for supported + options. Works for any key of :attr:`calibrations`. + + Returns + ------- + name : str + The generated name. + """ + return f"{self.cam.name}_{self.slm.name}_{calibration_type}" + + def save_calibration(self, calibration_type, path=".", name=None): + """ + to a file like ``"path/name_id.h5"``. + + Parameters + ---------- + calibration_type : str + The type of calibration to save. See :attr:`calibrations` for supported + options. Works for any key of :attr:`calibrations`. + path : str + Path to directory to save in. Default is current directory. + name : str OR None + Name of the save file. If ``None``, will use :meth:`name_calibration`. + + Returns + ------- + str + The file path that the calibration was saved to. + """ + if not calibration_type in self.calibrations: + raise ValueError(f"Could not find calibration '{calibration_type}' in calibrations.") + + if name is None: + name = self.name_calibration(calibration_type) + file_path = generate_path(path, name, extension="h5") + write_h5(file_path, {calibration_type : self.calibrations[calibration_type]}) + + return file_path + + def load_calibration(self, calibration_type, file_path=None): + """ + from a file. + + Parameters + ---------- + calibration_type : str + The type of calibration to save. See :attr:`calibrations` for supported + options. Works for any key of :attr:`calibrations`. + file_path : str OR None + Full path to the calibration file. If ``None``, will + search the current directory for a file with a name like + the one returned by :meth:`name_calibration`. + + Returns + ------- + str + The file path that the calibration was loaded from. + + Raises + ------ + FileNotFoundError + If a file is not found. + """ + if file_path is None: + path = os.path.abspath(".") + name = self.name_calibration(calibration_type) + file_path = latest_path(path, name, extension="h5") + if file_path is None: + raise FileNotFoundError( + "Unable to find a calibration file like\n{}" + "".format(os.path.join(path, name)) + ) + + self.calibrations[calibration_type] = read_h5(file_path)[calibration_type] + + return file_path + + def _get_calibration_metadata(self): + return { + "__version__" : __version__, + "__time__" : time.time(), + "__meta__" : { + "camera" : self.cam.name, # Future: pickle hardware. + "slm" : self.slm.name + } + } + ### Settle Time Calibration ### def settle_calibrate( - self, vector=(0.1, 0.1), basis="kxy", size=None, times=None, settle_time_s=1, plot=True + self, vector=(0.1, 0.1), basis="kxy", size=None, times=None, settle_time_s=1 ): """ - **(NotImplemented)** Approximates the :math:`1/e` settle time of the SLM. This is done by successively removing and applying a blaze to the SLM, - measuring the intensity at the first order spot versus time. + measuring the intensity at the first order spot versus time delay. Parameters ---------- vector : array_like - Point to TODO - basis : {"ij", "kxy"} - Basis of vector. This is the vector TODO + Point to measure settle time at via a simple blaze in the ``"kxy"`` basis. size : int - Size in pixels of the integration region. If ``None``, sets to sixteen - times the approximate size of a diffraction limited spot. - times : array_like + Size in pixels of the integration region in the ``"ij"`` basis. + If ``None``, sets to sixteen times the approximate size of a diffraction-limited spot. + times : array_like OR None OR int List of times to sweep over in search of the :math:`1/e` settle time. - settle_time_s : float - Time between measurements to allow the SLM to re-settle. - plot : bool - Whether to print debug plots. + If ``None``, defaults to 21 points over one second. + If an integer, defaults to that given number of points over one second. + settle_time_s : float OR None + Time between measurements to allow the SLM to re-settle. If ``None``, uses the + current default in the SLM. """ - if basis == "ij": - vector = self.ijcam_to_kxyslm(vector) - + # Parse vector. point = self.kxyslm_to_ijcam(vector) - blaze = toolbox.blaze(grid=self.slm, vector=vector) + # Parse size. + if size is None: + size = 16 * toolbox.convert_radius( + self.slm.spot_radius_kxy(), + to_units="ij", + slm=self + ) + size = int(size) + + # Parse times. + if times is None: + times = 21 + if np.isscalar(times): + times = np.linspace(0, 1, int(times), endpoint=True) + times = np.ravel(times) + + # Parse settle_time_s. + if settle_time_s is None: + settle_time_s = self.slm.settle_time_s + settle_time_s = float(settle_time_s) + results = [] + # Collect data for t in times: + # Reset the pattern and wait for it to settle self.slm.write(None, settle=False) time.sleep(settle_time_s) + # Turn on the pattern and wait for time t self.slm.write(blaze, settle=False) time.sleep(t) image = self.cam.get_image() + results.append(analysis.take(image, point, size, centered=True, integrate=True)) + + self.calibrations["settle"] = { + "times" : times, + "data" : results + } + self.calibrations["settle"].update(self._get_calibration_metadata()) - print(np.sum(image)) - print(np.max(image)) - # nonzero = image.ravel() - # nonzero = nonzero[nonzero != 0] - # print(nonzero) + self.process_settle_calibration(plot=False) - _, axs = plt.subplots(1, 2, figsize=(10, 4)) - axs[0].imshow(image) - axs[1].imshow( - np.squeeze(analysis.take(image, point, size, centered=True, integrate=False)) - ) - plt.show() + return self.calibrations["settle"] - results.append(analysis.take(image, point, size, centered=True, integrate=True)) + def process_settle_calibration(self, plot=True): + """ + Fits an exponential to the measured data to + approximate the :math:`1/e` settle time of the SLM. + + Parameters + ---------- + plot : bool + Whether to show a debug plot with the exponential fit. + + Returns + ------- + dict + The settle time and communication time measured. + """ + times = self.calibrations["settle"]["times"] + results = self.calibrations["settle"]["data"] if plot: plt.plot(times, np.squeeze(results), "k*") @@ -252,27 +399,105 @@ def settle_calibrate( plt.xlabel("Time [sec]") plt.show() - return results + times_clean = [] + results_clean = [] + + # TODO: Fit to step + exp instead. + + # Index of first non-null value + exp_start = -1 - def process_settle_calibration(self): - pass + # Number of consecutive non-null values for considering that the transient has started + to_check = 3 + + # Remove outliers + for i in range(len(results)): + # Check n consecutive elements, if they are non-null consider them valid from the previous one + if exp_start == -1 and i + to_check < len(results): + all_valid = True + for j in range(to_check): # Check the following numbers + if results[i+j] == 0: + all_valid = False + if all_valid == True: # Found n consecutive non-null elements + if i == 0: # Add the previous element to the results cleaned + exp_start = 0 + else: + exp_start = i-1 + results_clean.append(results[exp_start]) + times_clean.append(times[exp_start]) + results_clean.append(results[i]) + times_clean.append(times[i]) + else: + if results[i] != 0: # Remove zeros given by camera errors + results_clean.append(results[i]) + times_clean.append(times[i]) + + # Save communication delay, time tag before start of the transient + com_time = times[exp_start] + + # Function to interpolate + def exponential(x, a, b, c): + return c - a*np.exp(b*x) + + # Fit the date with the function + params, _ = optimize.curve_fit(exponential, times_clean, results_clean, maxfev = 10000) + a, b, c = params + + set_time = -1/b + + # Evaluate the fitting function in the interval + x_interp = np.linspace(min(times_clean), max(times_clean), 100) + y_interp = exponential(x_interp, a, b, c) + + if plot: + title = ( + f"Communication time: {int((1e3*com_time))} ms\n" + f"Settle time: {int((1e3*set_time))} ms" + ) + plt.plot(x_interp, y_interp, "--", linewidth=2, color='red', label='interpolation') + plt.plot(times, results, "k*", markersize=7, label='capta') + plt.xlabel("Time [sec]") + plt.ylabel("Signal [a.u.]") + plt.title(title) + plt.show() + + # Update dictionary with results. FUTURE: Return error bars? + processed = { + "settle_time" : set_time, + "communication_time" : com_time + } + self.calibrations["settle"].update(processed) + + return processed ### Pixel Crosstalk and Vpi Calibration ### def pixel_calibrate( self, - gray_levels=2, + levels=2, periods=2, - measure_orders=3, + orders=3, + window=None, + field_period=10, ): """ - Calibration for pixel properties such as crosstalk and . - Note this calibration is not at the level of individual pixels, but rather - aggregate statistics over all pixels. - - Phase crosstalk: https://doi.org/10.1364/OE.20.022334 - Vpi variation: https://doi.org/10.1364/OE.21.016086 - Upgraded model: https://doi.org/10.1364/OE.27.025046 + Physical SLMs do not produce perfectly sharp and discrete blocks of a desired + phase at each pixel. Rather, the realized phase might deviate from the desired + phase (error) and be blurred between pixels (crosstalk). + + We adopt a literature approach to calibrating both phenomena by `measuring the + system response of binary gratings `. + In the future, we intend to fit the measured data to `an upgraded asymmetric + model of phase crosstalk `_, and then + apply the model to beam propagation during holographic optimization. A better + understanding of the system error can lead to holograms that take this error + into account. + + Note that this algorithm does not operate at the level of individual pixels, but + rather on aggregate statistics over a region of pixels. + Right now, this calibration is done for one region (which defaults to the full + SLM). In the future, we might want to calibrate many regions across the SLM to + measure `spatially varying phase response `_ Note ~~~~ @@ -280,39 +505,63 @@ def pixel_calibrate( Caution ~~~~~~~ - Data is taken without wavefront calibration applied. If the SLM without - calibration produces too defocussed of a spot, - then this measurement may not be ideal. + Data must be acquired without wavefront calibration applied. + If the uncalibrated SLM produces too defocussed of a spot, + then this measurement may not be ideal. On the flip side, a + too-focussed spot might increase error by integration over fewer camera pixels. Parameters ---------- - gray_levels : int OR array_like of int - TODO + levels : int OR array_like of int + Which bitlevels to test, out of the :math:`2^B` levels available for a + :math:`B`-bit SLM. Note that runtime scales with :math:`\mathcal{O}(B^2)`. + If an integer is passed, the integer is rounded to the next largest power of + two, and this number of bitlevels are sampled. periods : int OR array_like of int - TODO - measure_orders : int OR array_like of int - TODO - plot : bool - TODO + Periods (in pixels) of the binary gratings that we will apply. + Must be even integers. + orders : int OR array_like of int + Orders (..., -1st, 0th, 1st, ...) of the binary gratings to measure data at. + If scalar is provided, measures orders between -nth and nth order, inclusive. + window + If not ``None``, the pixel calibration is only done over the region of the SLM + defined by ``window``. + Passed to :meth:`~slmsuite.holography.toolbox.window_slice()`. + See :meth:`~slmsuite.holography.toolbox.window_slice()` for various options. + field_period : int + If ``window`` is not ``None``, then the field is deflected away in an + orthogonal direction with a grating of the given period. """ - # Parse gray_levels by forcing range and datatype. - if np.isscalar(gray_levels): - gray_levels = np.todo() - gray_levels = np.mod(gray_levels, self.slm.bitresolution).astype(self.slm.dtype) - N = len(gray_levels) + # Parse levels by forcing range and datatype. + if np.isscalar(levels): + if levels < 1: + levels = 1 + levels = 2 ** (np.ceil(np.log2(levels))) - # Parse periods by forcing even integer. - if np.isscalar(measure_orders): - measure_orders = int(measure_orders) - measure_orders = np.arange(-measure_orders, measure_orders+1) + if levels > self.slm.bitresolution: + warnings.warn("Requested more levels than available. Rounding down.") + levels = self.slm.bitresolution - measure_orders = measure_orders.astype(int) - M = len(measure_orders) + levels = np.arange(levels) * (self.slm.bitresolution / levels) + levels = np.mod(levels, self.slm.bitresolution).astype(self.slm.dtype) + N = len(levels) - # Parse periods by forcing even integer. + # Parse periods by forcing integer. periods = periods.astype(int) + periods = 2 * (periods // 2) P = len(periods) + if np.any(periods <= 0): + raise ValueError("period should not be negative.") + + # Parse orders by forcing integer. + if np.isscalar(orders): + orders = int(orders) + orders = np.arange(-orders, orders+1) + orders = orders.astype(int) + M = len(orders) + + # Figure out our shape. shape = (2, P, N, N, M) data = np.zeros(shape) @@ -326,18 +575,40 @@ def pixel_calibrate( slm=self ) + # Make the y-pointing field vector, then the x-pointing field vector. + field_freq = np.zeros((2, 2)) + field_freq[0, 0] = field_freq[1, 1] = np.reciprocal(field_period.astype(float)) + field_kxy = toolbox.convert_vector( + field_freq, + from_units="freq", + to_units="norm", + slm=self + ) + field_hi = self.slm.dtype(self.slm.bitresolution / 2) + field_lo = self.slm.dtype(0) + + field_ij = toolbox.convert_vector( + field_freq, + from_units="freq", + to_units="ij", + slm=self + ) + # Figure out where the orders will appear on the camera. vectors_ij = self.kxyslm_to_ijcam(vectors_kxy) center = self.kxyslm_to_ijcam((0,0)) dorder = vectors_ij - center + dfield = field_ij - center order_ij = [] for i in range(2*P): - order_ij.append(center + measure_orders * dorder[:, [i]]) + order_ij.append(center + orders * dorder[:, [i]]) - distances = np.linalg.norm(dorder, axis=0) - integration_size = np.ceil(np.min(distances)) + integration_size = int(np.ceil(np.min([ + np.min(np.max(dorder, axis=1)), + np.min(np.max(dfield, axis=1)) + ]))) # Warn the user if any order is outside the field of view. if False: @@ -349,12 +620,34 @@ def pixel_calibrate( for j in range(P): # Period for k in range(N): # Upper triangular gray level selection. for l in range(k, N): # Periodic normalization when equal. - phase = toolbox.phase.binary( - self.slm, - vectors_kxy[:, prange[j]], - gray_levels[k], - gray_levels[l] - ) + if window is None: + phase = toolbox.phase.binary( + self.slm, + vector=vectors_kxy[:, prange[j]], + a=levels[k], + b=levels[l] + ) + else: + # In windowed mode, blaze the field away from the 0th order, + # in the direction perpendicular to the target. + phase = toolbox.phase.binary( + grid=self.slm, + vector=field_kxy[:, i], + a=field_hi, + b=field_lo + ) + toolbox.imprint( + phase, + window=window, + function=toolbox.phase.binary, + grid=self.slm, + vector=vectors_kxy[:, prange[j]], + a=levels[k], + b=levels[l] + ) + + # We're writing integers, so this goes directly to the SLM, + # bypassing phase2gray. self.slm.write(phase, phase_calibrate=False, settle=True) data[i,j,k,l,:] = data[i,j,l,k,:] = analysis.take( @@ -364,8 +657,36 @@ def pixel_calibrate( integrate=True, ).astype(float) + # Assemble the return dictionary. + self.calibrations["pixel"] = { + "levels" : levels, + "periods" : periods, + "orders" : orders, + "data": data + } + self.calibrations["pixel"].update(self._get_calibration_metadata()) + + # Process by default because we currently don't have any arguments. + self.process_pixel_calibration() + + return self.calibrations["pixel"] + def process_pixel_calibration(self): - pass + """ + Currently, this method only displays debug plots of the measurements. + In the future, the measurements will be fit in a way that can be applied to + propagation. + """ + periods = self.calibrations["pixel"]["periods"] + orders = self.calibrations["pixel"]["orders"] + data = self.calibrations["pixel"]["data"] + + for i, direction in enumerate(["x", "y"]): + for j, period in enumerate(periods): + for o, order in enumerate(orders): + plt.imshow(data[i,j,:,:,o]) + plt.title(f"{period}-pixel, ${direction}$ grating; measuring order {order}") + plt.show() ### Fourier Calibration ### @@ -382,7 +703,6 @@ def fourier_calibrate( """ Project and fit a SLM computational Fourier space ``"knm"`` grid onto camera pixel space ``"ij"`` for affine fitting. - Sets :attr:`~slmsuite.hardware.cameraslms.FourierSLM.fourier_calibration`. An array produced by :meth:`~slmsuite.holography.algorithms.SpotHologram.make_rectangular_array()` is projected for analysis by @@ -427,7 +747,7 @@ def fourier_calibrate( Returns ------- dict - :attr:`~slmsuite.hardware.cameraslms.FourierSLM.fourier_calibration` + :attr:`~slmsuite.hardware.cameraslms.FourierSLM.calibrations["fourier"]` """ # Parse variables if isinstance(array_shape, REAL_TYPES): @@ -498,83 +818,18 @@ def fourier_calibrate( [M[1, 0] * scaling[0], M[1, 1] * scaling[1]], ]) - self.fourier_calibration = {"__version__" : __version__, "M": M, "b": b, "a": a} + self.calibrations["fourier"] = { + "__version__" : __version__, + "__time__" : time.time(), + "M": M, + "b": b, + "a": a + } - return self.fourier_calibration + return self.calibrations["fourier"] ### Fourier Calibration Helpers ### - def name_fourier_calibration(self): - """ - Creates ``"--fourier-calibration"``. - - Returns - ------- - name : str - The generated name. - """ - return "{}-{}-fourier-calibration".format(self.cam.name, self.slm.name) - - def save_fourier_calibration(self, path=".", name=None): - """ - Saves :attr:`~slmsuite.hardware.cameraslms.FourierSLM.fourier_calibration` - to a file like ``"path/name_id.h5"``. - - Parameters - ---------- - path : str - Path to directory to save in. Default is current directory. - name : str OR None - Name of the save file. If ``None``, will use :meth:`name_fourier_calibration`. - - Returns - ------- - str - The file path that the Fourier calibration was saved to. - """ - if name is None: - name = self.name_fourier_calibration() - file_path = generate_path(path, name, extension="h5") - write_h5(file_path, self.fourier_calibration) - - return file_path - - def load_fourier_calibration(self, file_path=None): - """ - Loads :attr:`~slmsuite.hardware.cameraslms.FourierSLM.fourier_calibration` - from a file. - - Parameters - ---------- - file_path : str OR None - Full path to the Fourier calibration file. If ``None``, will - search the current directory for a file with a name like - the one returned by :meth:`name_fourier_calibration`. - - Returns - ------- - str - The file path that the Fourier calibration was loaded from. - - Raises - ------ - FileNotFoundError - If a file is not found. - """ - if file_path is None: - path = os.path.abspath(".") - name = self.name_fourier_calibration() - file_path = latest_path(path, name, extension="h5") - if file_path is None: - raise FileNotFoundError( - "Unable to find a Fourier calibration file like\n{}" - "".format(os.path.join(path, name)) - ) - - self.fourier_calibration = read_h5(file_path) - - return file_path - def project_fourier_grid(self, array_shape=10, array_pitch=10, array_center=None, **kwargs): """ Projects a Fourier space grid ``"knm"`` onto pixel space ``"ij"``. @@ -641,7 +896,12 @@ def set_fourier_calibration_analytic(self, M, b): a = format_2vectors([0,0]) b = format_2vectors(b) - self.fourier_calibration = {"__version__" : __version__, "M": M, "b": b, "a": a} + self.calibrations["fourier"] = { + "M": M, + "b": b, + "a": a + } + self.calibrations["fourier"].update(self._get_calibration_metadata()) def build_fourier_calibration( self, @@ -653,7 +913,7 @@ def build_fourier_calibration( ): """ Meta: This docstring will be overwritten by ``SimulatedCamera.build_affine``'s - at the end of this file. + after this class. """ if offset is None: offset = np.flip(self.cam.shape) / 2 @@ -695,7 +955,6 @@ def kxyslm_to_ijcam(self, kxy): .. math:: \vec{y} = M \cdot (\vec{x} - \vec{a}) + \vec{b} where :math:`M`, :math:`\vec{b}`, and :math:`\vec{a}` are stored in - :attr:`~slmsuite.hardware.cameraslms.FourierSLM.fourier_calibration`. If the vectors are three-dimensional, the third depth dimension is treated according to: @@ -725,16 +984,16 @@ def kxyslm_to_ijcam(self, kxy): RuntimeError If the fourier plane calibration does not exist. """ - if self.fourier_calibration is None: + if self.calibrations["fourier"] is None: raise RuntimeError("Fourier calibration must exist to be used.") kxy = format_vectors(kxy, handle_dimension="pass") # Apply the xy transformation. ij = np.matmul( - self.fourier_calibration["M"], - kxy[:2, :] - self.fourier_calibration["a"] - ) + self.fourier_calibration["b"] + self.calibrations["fourier"]["M"], + kxy[:2, :] - self.calibrations["fourier"]["a"] + ) + self.calibrations["fourier"]["b"] # Handle z if needed. if kxy.shape[0] == 3: @@ -750,7 +1009,7 @@ def ijcam_to_kxyslm(self, ij): .. math:: \vec{x} = M^{-1} \cdot (\vec{y} - \vec{b}) + \vec{a} where :math:`M`, :math:`\vec{b}`, and :math:`\vec{a}` are stored in - :attr:`~slmsuite.hardware.cameraslms.FourierSLM.fourier_calibration`. + :attr:`~slmsuite.hardware.cameraslms.FourierSLM.calibrations["fourier"]`. Important ~~~~~~~~~ @@ -787,16 +1046,16 @@ def ijcam_to_kxyslm(self, ij): RuntimeError If the fourier plane calibration does not exist. """ - if self.fourier_calibration is None: + if self.calibrations["fourier"] is None: raise RuntimeError("Fourier calibration must exist to be used.") ij = format_vectors(ij, handle_dimension="pass") # Apply the xy transformation. kxy = np.matmul( - np.linalg.inv(self.fourier_calibration["M"]), - ij[:2, :] - self.fourier_calibration["b"] - ) + self.fourier_calibration["a"] + np.linalg.inv(self.calibrations["fourier"]["M"]), + ij[:2, :] - self.calibrations["fourier"]["b"] + ) + self.calibrations["fourier"]["a"] # Handle z if needed. if ij.shape[0] == 3: @@ -814,7 +1073,7 @@ def get_farfield_spot_size(self, slm_size=None, basis="kxy"): is accomplished using the calibration produced by :meth:`~slmsuite.hardware.cameraslms.FourierSLM.fourier_calibrate()` and stored in - :attr:`~slmsuite.hardware.cameraslms.FourierSLM.fourier_calibration`. + :attr:`~slmsuite.hardware.cameraslms.FourierSLM.calibrations["fourier"]`. Parameters ---------- @@ -846,7 +1105,7 @@ def get_farfield_spot_size(self, slm_size=None, basis="kxy"): if basis == "kxy": return (1 / slm_size[0], 1 / slm_size[1]) elif basis == "ij": - M = self.fourier_calibration["M"] + M = self.calibrations["fourier"]["M"] # Compensate for spot rotation s.t. spot size is along camera axes size_kxy = np.linalg.inv(M / np.sqrt(np.linalg.det(M))) @ np.array( (1 / slm_size[0], 1 / slm_size[1]) @@ -885,11 +1144,11 @@ def get_effective_focal_length(self, units="norm"): f_eff : float Effective focal length. """ - if self.fourier_calibration is None: + if self.calibrations["fourier"] is None: raise RuntimeError("Fourier calibration must exist to be used.") # Gather f_eff in pix/rad. - f_eff = np.sqrt(np.abs(np.linalg.det(self.fourier_calibration["M"]))) + f_eff = np.sqrt(np.abs(np.linalg.det(self.calibrations["fourier"]["M"]))) # Gather other conversions. if units != "ij" and self.cam.pitch_um is None: @@ -910,6 +1169,138 @@ def get_effective_focal_length(self, units="norm"): ### Wavefront Calibration ### def wavefront_calibrate( + self, + *args, + method=None, + **kwargs, + ): + if method is None: + method = "superpixel" + + if method == "superpixel": + return self.wavefront_calibrate_superpixel(*args, **kwargs) + elif method == "zernike": + return self.wavefront_calibrate_zernike(*args, **kwargs) + else: + raise ValueError(f"Wavefront calibration method {method} not recognized.") + + ### Zernike Wavefront Calibration ### + + def wavefront_calibrate_zernike( + self, + calibration_points=None, + zernike_indices=3, + wavefronts=1, + phase_steps=11, + callback=None, + metric=None, + plot=0, + ): + """ + Perform wavefront calibration. + + Parameters + ---------- + calibration_points : (float, float) OR numpy.ndarray OR None + Position(s) in the camera domain where interference occurs. + For multiple positions, this must be of shape ``(2, N)``. + This is naturally in the ``"ij"`` basis. + If ``None``, densely fills the camera field of view with calibration points. + zernike_indices : int OR list of int + TODO + wavefronts : float OR list of float + TODO + plot : int or bool + Whether to provide visual feedback, options are: + + - ``-1`` : No plots or tqdm prints. + - ``0``, ``False`` : No plots, but tqdm prints. + - ``1``, ``True`` : Plots on fits and essentials. + - ``2`` : Plots on everything. + - ``3`` : ``test_index`` not ``None`` only: returns image frames + to make a movie from the phase measurement (not for general use). + + Returns + ------- + dict + The contents of + :attr:`~slmsuite.hardware.cameraslms.FourierSLM.calibrations["wavefront"]`. + + Raises + ------ + RuntimeError + If the Fourier plane calibration does not exist. + ValueError + If various points are out of range. + """ + + + # def get_default_callback(cam, points, window_size, metric=None): + + def sweep_me(slm, sweep, perturbation, pattern, callback): + result = None + sweep = np.ravel(sweep) + N = len(sweep) + M = None + + for i, x in enumerate(sweep): + phase = pattern + x * perturbation + slm.write(phase) + this_result = np.array(callback()) + + if result is None: + M = len(this_result) + result = np.full((N, M), np.nan, dtype=this_result.dtype) + + if len(this_result) != M: + raise RuntimeError() + else: + result[i, :] = this_result + + return result + + if np.isscalar(zernike_indices): + # Default to the lowest indices + # + 3 to avoid the piston, x tilt, y tilt + zernike_indices = np.arange(zernike_indices) + 3 + + + window_size = None + + if callback is None: + def default_callback(): + self.cam.flush() + img = self.cam.get_image() + + images = analysis.take(img, calibration_points, window_size) + + if metric is None: + return FourierSLM._wavefront_calibrate_zernike_default_metric(images) + else: + return metric(images) + + callback = default_callback + + sweep = np.linspace(-wavefronts, wavefronts, phase_steps, endpoint=True) + + for i in zernike_indices: + perturbation = toolbox.phase.zernike(self.slm, i) + + result = sweep_me(self.slm, sweep, perturbation, pattern, callback) + + # analyze + + @staticmethod + def _wavefront_calibrate_zernike_default_metric(images): + """ + Stack of images + """ + variances = analysis.image_variances(images) + return analysis.image_areas(variances) + + ### Superpixel Wavefront Calibration ### + + def wavefront_calibrate_superpixel( self, calibration_points=None, superpixel_size=50, @@ -935,7 +1326,7 @@ def wavefront_calibrate( from each point, the less ideal the correction is. Correction at many points over the plane permits a better understanding of the aberration and greater possibility of compensation. - Sets :attr:`~slmsuite.hardware.cameraslms.FourierSLM.wavefront_calibration_raw`. + Sets :attr:`~slmsuite.hardware.cameraslms.FourierSLM.calibrations["wavefront"]`. Run :meth:`~slmsuite.hardware.cameraslms.FourierSLM.process_wavefront_calibration` after to produce the usable calibration which is written to the SLM. This procedure measures the wavefront phase and amplitude. @@ -1021,7 +1412,7 @@ def wavefront_calibrate( calibration. This is useful to determine the quality of a previous calibration, as a new calibration should yield zero phase correction needed if the previous was perfect. The old calibration will be stored in the - :attr:`wavefront_calibration_raw` under ``"previous_phase_correction"``, + :attr:`calibrations` under ``"previous_phase_correction"``, so keep in mind that this (uncompressed) image will take up significantly more space. measure_background : bool @@ -1044,7 +1435,7 @@ def wavefront_calibrate( ------- dict The contents of - :attr:`~slmsuite.hardware.cameraslms.FourierSLM.wavefront_calibration_raw`. + :attr:`~slmsuite.hardware.cameraslms.FourierSLM.calibrations["wavefront"]`. Raises ------ @@ -1139,7 +1530,7 @@ def index2image(index): field_point = np.rint(format_2vectors(field_point)).astype(int) # Use the Fourier calibration to help find points/sizes in the imaging plane. - if self.fourier_calibration is None: + if self.calibrations["fourier"] is None: raise RuntimeError("Fourier calibration must be done before wavefront calibration.") calibration_blazes = self.ijcam_to_kxyslm(calibration_points) reference_blazes = calibration_blazes.copy() @@ -1318,6 +1709,7 @@ def index2image(index): # Build the calibration dict. calibration_dict = { "__version__" : __version__, + "__time__" : time.time(), "calibration_points" : calibration_points, "superpixel_size" : superpixel_size, "slm_supershape" : slm_supershape, @@ -1511,7 +1903,7 @@ def fit_phase_image(img, dsuperpixel, plot_fits=True): xyr = [l.ravel() for l in xy] # Process dsuperpixel by rotating it according to the Fourier calibration. - M = self.fourier_calibration["M"] + M = self.calibrations["fourier"]["M"] M_norm = M / np.sqrt(np.abs(np.linalg.det(M))) dsuperpixel = np.squeeze(np.matmul(M_norm, format_2vectors(dsuperpixel))) @@ -1982,89 +2374,13 @@ def measure(schedule, plot=False): calibration_dict[key][i, coords[1, i], coords[0, i]] = result - self.wavefront_calibration_raw = calibration_dict + self.calibrations["wavefront"] = calibration_dict + self.calibrations["wavefront"].update(self._get_calibration_metadata()) return calibration_dict ### Wavefront Calibration Helpers ### - def name_wavefront_calibration(self): - """ - Creates ``"--wavefront-calibration"``. - - Returns - ------- - name : str - The generated name. - """ - return "{}-{}-wavefront-calibration".format(self.cam.name, self.slm.name) - - def save_wavefront_calibration(self, path=".", name=None): - """ - Saves :attr:`~slmsuite.hardware.cameraslms.FourierSLM.wavefront_calibration_raw` - to a file like ``"path/name_id.h5"``. - - Parameters - ---------- - path : str - Path to the save location. Default is current directory. - name : str or None - Name of the save file. If ``None``, will use :meth:`name_wavefront_calibration`. - - Returns - ------- - str - The file path that the fourier calibration was saved to. - """ - if name is None: - name = self.name_wavefront_calibration() - file_path = generate_path(path, name, extension="h5") - write_h5(file_path, self.wavefront_calibration_raw) - - return file_path - - def load_wavefront_calibration(self, file_path=None, process=True, **kwargs): - """ - Loads :attr:`~slmsuite.hardware.cameraslms.FourierSLM.wavefront_calibration_raw` - from a file. - - Parameters - ---------- - file_path : str or None - Full path to the wavefront calibration file. If ``None``, will - search the current directory for a file with a name like - the one returned by - :meth:`~slmsuite.hardware.cameraslms.FourierSLM.name_wavefront_calibration`. - process : bool - Whether to immediately process the wavefront calibration. - See - :meth:`~slmsuite.hardware.cameraslms.FourierSLM.process_wavefront_calibration`. - **kwargs - Passed to :meth:`~slmsuite.hardware.cameraslms.FourierSLM.process_wavefront_calibration`, - if ``process`` is true. - - Returns - ------- - str - The file path that the wavefront calibration was loaded from. - """ - if file_path is None: - path = os.path.abspath(".") - name = self.name_wavefront_calibration() - file_path = latest_path(path, name, extension="h5") - if file_path is None: - raise FileNotFoundError( - "Unable to find a Fourier calibration file like\n{}" - "".format(os.path.join(path, name)) - ) - - self.wavefront_calibration_raw = read_h5(file_path) - - if process: - self.process_wavefront_calibration(**kwargs) - - return file_path - def get_wavefront_calibration_points( self, pitch, @@ -2219,7 +2535,7 @@ def get_wavefront_calibration_window(self, superpixel_size): def process_wavefront_calibration(self, index=0, smooth=True, r2_threshold=0.9, apply=True, plot=False): """ - Processes :attr:`~slmsuite.hardware.cameraslms.FourierSLM.wavefront_calibration_raw` + Processes :attr:`~slmsuite.hardware.cameraslms.FourierSLM.calibrations["wavefront"]` into the desired phase correction and amplitude measurement. Applies these parameters to the respective variables in the SLM if ``apply`` is ``True``. @@ -2249,7 +2565,7 @@ def process_wavefront_calibration(self, index=0, smooth=True, r2_threshold=0.9, The updated source dictionary containing the processed source amplitude and phase. """ # Step 0: Initialize helper variables and functions. - data = self.wavefront_calibration_raw + data = self.calibrations["wavefront"] if len(data) == 0: raise RuntimeError("No raw wavefront data to process. Either load data or calibrate.") @@ -2299,7 +2615,7 @@ def _process_wavefront_calibration_r001(self, smooth=True, r2_threshold=0.9, app The updated source dictionary containing the processed source amplitude and phase. """ # Step 0: Initialize helper variables and functions. - data = self.wavefront_calibration_raw + data = self.calibrations["wavefront"] if len(data) == 0: raise RuntimeError("No raw wavefront data to process. Either load data or calibrate.") @@ -2548,7 +2864,7 @@ def plot_wavefront_calibration_raw(self, index=None, r2_threshold=0, phase_detai """TODO""" plt.figure(figsize=(16, 8)) - data = self.wavefront_calibration_raw + data = self.calibrations["wavefront"] if index is None: coords = data["calibration_points"] @@ -2639,4 +2955,14 @@ def plot_wavefront_calibration_raw(self, index=None, r2_threshold=0, phase_detai plt.show() -FourierSLM.build_fourier_calibration.__doc__ = SimulatedCamera.build_affine.__doc__ \ No newline at end of file +FourierSLM.build_fourier_calibration.__doc__ = SimulatedCamera.build_affine.__doc__ + + +class _SimulatedFourierSLM(FourierSLM): + def __init__(calibrations): + """ + Parameters + ---------- + calibrations : list of str + TODO + """ diff --git a/slmsuite/hardware/slms/slm.py b/slmsuite/hardware/slms/slm.py index 0980c78..77838e9 100644 --- a/slmsuite/hardware/slms/slm.py +++ b/slmsuite/hardware/slms/slm.py @@ -755,7 +755,7 @@ def fit_source_amplitude(self, method="fit", extent_threshold=.1, force=True): return center_grid = np.array( - [np.argmin(self.grid[0][0,:]), np.argmin(self.grid[1][:,0])] + [np.argmin(np.abs(self.grid[0][0,:])), np.argmin(np.abs(self.grid[1][:,0]))] ) if not "amplitude" in self.source: @@ -780,21 +780,22 @@ def fit_source_amplitude(self, method="fit", extent_threshold=.1, force=True): raise RuntimeError("extent_threshold cannot exceed 1 (100%). Use a small value.") if method == "fit": - result = analysis.image_fit(amp, plot=True) - + result = analysis.image_fit(amp, plot=False) std = np.array([result[0,5], result[0,6]]) center = np.array([result[0,1], result[0,2]]) - center += np.flip(self.shape)/2 elif method == "moments": center = analysis.image_positions(amp) - std = np.sqrt(analysis.image_variances(amp, centers=center)[:2]) - # std_norm = self.pitch * np.squeeze(std_pix) - # std = np.mean(std_norm) center = np.squeeze(center) - center += np.flip(self.shape)/2 + + print(center) + + center += np.flip(self.shape)/2 + + + print(center) self.source["amplitude_center_pix"] = center self.source["amplitude_radius"] = np.mean(self.pitch * np.squeeze(std)) @@ -802,8 +803,17 @@ def fit_source_amplitude(self, method="fit", extent_threshold=.1, force=True): # Handle centering. dcenter = center_grid - center - self.grid[0] -= dcenter[0] * self.pitch[0] - self.grid[1] -= dcenter[1] * self.pitch[1] + print(dcenter) + print(self.pitch) + + self.grid[0] += dcenter[0] * self.pitch[0] + self.grid[1] += dcenter[1] * self.pitch[1] + + + center_grid = np.array( + [np.argmin(self.grid[0][0,:]), np.argmin(self.grid[1][:,0])] + ) + print(center_grid) extent_mask = amp > (extent_threshold * np.amax(amp)) @@ -817,7 +827,7 @@ def fit_source_amplitude(self, method="fit", extent_threshold=.1, force=True): def get_source_radius(self): """ - Extracts the source radius for + Extracts the source radius in normalized units for functions like :meth:`~slmsuite.holography.toolbox.phase.laguerre_gaussian()` from the scalars computed in :meth:`~slmsuite.hardware.slms.slm.SLM.fit_source_amplitude()`. @@ -872,7 +882,7 @@ def plot_source(self, sim=False): elif not sim and not np.all([k in self.source for k in ("amplitude", "phase")]): raise RuntimeError( "'amplitude' or 'phase' keywords missing from slm.source! Run " - ".wavefront_calibration() or .set_source_analytic() to measure source profile." + ".wavefront_calibrate() or .set_source_analytic() to measure source profile." ) plot_r2 = not sim and "r2" in self.source diff --git a/slmsuite/holography/algorithms/__init__.py b/slmsuite/holography/algorithms/__init__.py index 010f892..f97351d 100644 --- a/slmsuite/holography/algorithms/__init__.py +++ b/slmsuite/holography/algorithms/__init__.py @@ -16,56 +16,11 @@ Tip ~~~ This module makes use of the GPU-accelerated computing library :mod:`cupy` -(`GitHub `_) +(`GitHub `_). If :mod:`cupy` is not supported, then :mod:`numpy` is used as a fallback, though CPU alone is significantly slower. Using :mod:`cupy` is highly encouraged. - -Important ---------- -:mod:`slmsuite` follows the ``shape = (h, w)`` and ``vector = (x, y)`` formalism adopted by -the :mod:`numpy` ecosystem. :mod:`numpy`, :mod:`scipy`, :mod:`matplotlib`, etc generally follow this -formalism. The ``shape`` and indexing of an array or image always uses the inverted ``(h, w)`` form, -but other functions such as ``numpy.meshgrid(x, y)`` (default), ``scipy.odr.Data(x, y)``, or -``matplotlib.pyplot.scatter(x, y)`` use the standard cartesian ``(x, y)`` form that is more familiar -to users. This is not ideal and causes confusion, but this is the formalism generally -adopted by the community. - -Important -~~~~~~~~~ -This package uses a number of bases or coordinate spaces. Some coordinate spaces are -directly used by the user (most often the camera basis ``"ij"`` used for feedback). -Other bases are less often used directly, but are important to how holograms are -optimized under the hood (esp. ``"knm"``, the coordinate space of optimization when -using discrete Fourier transforms). - -.. list-table:: Bases used in :mod:`slmsuite`. - :widths: 20 80 - :header-rows: 1 - - * - Basis - - Meaning - * - ``"ij"`` - - Pixel basis of the camera. Centered at ``(i, j) = (cam.shape[1]/2, - cam.shape[0]/2)``. Is in the image space of the camera. - * - ``"kxy"`` - - Normalized (floating point) basis of the SLM's :math:`k`-space in normalized units. - Centered at ``(kx, ky) = (0, 0)``. This basis is what the SLM projects in angular - space (which maps to the camera's image space via the Fourier transform implemented by free - space and solidified by a lens). - * - ``"knm"`` - - Pixel basis of the SLM's computational :math:`k`-space. Centered at ``(kn, km) = - (shape[1]/2, shape[0]/2)``. - ``"knm"`` is a discrete version of the continuous ``"kxy"``. This is - important because holograms need to be stored in computer memory, a discrete - medium with pixels, rather than being purely continuous. For instance, in - :class:`SpotHologram`, spots targeting specific continuous angles are rounded to - the nearest discrete pixels of ``"knm"`` space in practice. - Then, this ``"knm"`` space image is handled as a - standard image/array, and operations such as the discrete Fourier transform - (instrumental for numerical hologram optimization) can be applied. - -See the first tip in :class:`Hologram` to learn more about ``"kxy"`` and ``"knm"`` -space. +The only algorithm that has no CPU fallback is :class:`CompressedSpotHologram` +in a Zernike basis beyond 2D or 3D spots. Other cases can use the CPU (however slowly). Note ~~~~ @@ -73,10 +28,11 @@ to enhance clarity and reduce file length. - ``_header.py`` : The common imports for all the files. -- ``_hologram.py`` : The core file. Contains the actual algorithms. -- ``_stats.py`` : Statistics and plotting common to all Holograms. -- ``_feedback.py`` : Infrastructure for image feedback. -- ``_spots.py`` : Infrastructure for spot-specific holography. +- ``_stats.py`` : Statistics and plotting common to all hologram classes. +- ``_hologram.py`` : The core file. Contains the actual algorithms (:class:`Hologram`). +- ``_feedback.py`` : Infrastructure for image feedback (:class:`FeedbackHologram`). +- ``_spots.py`` : Infrastructure for spot-specific holography + (:class:`SpotHologram`, :class:`CompressedSpotHologram`). """ from slmsuite.holography.algorithms._header import * diff --git a/slmsuite/holography/algorithms/_feedback.py b/slmsuite/holography/algorithms/_feedback.py index fe3c063..99e13da 100644 --- a/slmsuite/holography/algorithms/_feedback.py +++ b/slmsuite/holography/algorithms/_feedback.py @@ -10,14 +10,13 @@ class FeedbackHologram(Hologram): Attributes ---------- - cameraslm : slmsuite.hardware.cameraslms.FourierSLM OR None + cameraslm : slmsuite.hardware.cameraslms.FourierSLM OR None OR (int, int) A hologram with experimental feedback needs access to an SLM and camera. - If None, no feedback is applied (mostly defaults to :class:`Hologram`). - cam_shape : (int, int) - Shape of the camera in the meaning of :meth:`numpy.shape()`. - cam_points : numpy.ndarray - Array containing points corresponding to the corners of the camera in the SLM's k-space. - First point is repeated at the end for easy plotting. + If ``None``, no feedback is applied (mostly defaults to :class:`Hologram`). + _cam_points : numpy.ndarray + Array containing points corresponding to the corners of the camera in the SLM's + k-space. At the moment, this is only used for plotting. + First point is repeated at the end. target_ij : array_like OR None Amplitude target in the ``"ij"`` (camera) basis. Of same ``shape`` as the camera in :attr:`cameraslm`. Counterpart to :attr:`target` which is in the ``"knm"`` @@ -38,16 +37,16 @@ def __init__(self, shape, target_ij=None, cameraslm=None, **kwargs): shape : (int, int) Computational shape of the SLM in :mod:`numpy` `(h, w)` form. See :meth:`.Hologram.__init__()`. target_ij : array_like OR None - See :attr:`target_ij`. Should only be ``None`` if the :attr:`target` - will be generated by other means (see :class:`SpotHologram`), so the - user should generally provide an array. + See :attr:`target_ij`. + Is ``None`` only if the :attr:`target` will be generated by other means + (used by subclass :class:`SpotHologram`). This should not generally be used. cameraslm : slmsuite.hardware.cameraslms.FourierSLM OR slmsuite.hardware.slms.SLM OR None Provides access to experimental feedback. - If an :class:`slmsuite.hardware.slms.SLM` is passed, this is set to `None`, + If an :class:`~slmsuite.hardware.slms.SLM` is passed, the attribute is set to ``None``, but the information contained in the SLM is passed to the superclass :class:`.Hologram`. See :attr:`cameraslm`. **kwargs - See :meth:`Hologram.__init__`. + Passed to :meth:`Hologram.__init__`. """ # Use the Hologram constructor to initialize self.target with proper shape, # pass other arguments (esp. slm_shape). @@ -85,7 +84,7 @@ def __init__(self, shape, target_ij=None, cameraslm=None, **kwargs): else: self.target_ij = target_ij.astype(self.dtype) - if self.cameraslm is not None and self.cameraslm.fourier_calibration is not None: + if self.cameraslm is not None and "fourier" in self.cameraslm.calibrations: # Generate a list of the corners of the camera, for plotting. cam_shape = self.cameraslm.cam.shape @@ -96,18 +95,16 @@ def __init__(self, shape, target_ij=None, cameraslm=None, **kwargs): points_ij = toolbox.format_2vectors(np.vstack((ll, lr, ur, ul, ll)).T) points_kxy = self.cameraslm.ijcam_to_kxyslm(points_ij) - self.cam_points = toolbox.convert_vector( + self._cam_points = toolbox.convert_vector( points_kxy, "kxy", "knm", slm=self.cameraslm.slm, shape=self.shape ) - self.cam_shape = cam_shape # Transform the target, if it is provided. if target_ij is not None: self.update_target(target_ij, reset_weights=True) else: - self.cam_points = None - self.cam_shape = None + self._cam_points = None # Image transformation helper function. def ijcam_to_knmslm(self, img, out=None, blur_ij=None, order=3): @@ -139,23 +136,24 @@ def ijcam_to_knmslm(self, img, out=None, blur_ij=None, order=3): numpy.ndarray OR cupy.ndarray Image transformed into ``"knm"`` space. """ - assert self.cameraslm is not None - assert self.cameraslm.fourier_calibration is not None - - # First transformation. - conversion = toolbox.convert_vector( - (1, 1), "knm", "kxy", slm=self.cameraslm.slm, shape=self.shape - ) - toolbox.convert_vector( - (0, 0), "knm", "kxy", slm=self.cameraslm.slm, shape=self.shape + if self.cameraslm is None: + raise RuntimeError("Cannot use ijcam_to_knmslm without the calibrations in a cameraslm.") + if not "fourier" in self.cameraslm.calibrations: + raise RuntimeError("ijcam_to_knmslm requires a Fourier calibration.") + + # First transformation. FUTURE: make convert_basis output a matrix? + conversion = ( + toolbox.convert_vector((1, 1), "knm", "kxy", slm=self.cameraslm.slm, shape=self.shape) - + toolbox.convert_vector((0, 0), "knm", "kxy", slm=self.cameraslm.slm, shape=self.shape) ) M1 = np.diag(np.squeeze(conversion)) b1 = np.matmul(M1, -toolbox.format_2vectors(np.flip(np.squeeze(self.shape)) / 2)) # Second transformation. - M2 = self.cameraslm.fourier_calibration["M"] - b2 = self.cameraslm.fourier_calibration["b"] - np.matmul( - M2, self.cameraslm.fourier_calibration["a"] - ) + M2 = self.cameraslm.calibrations["fourier"]["M"] + b2 = self.cameraslm.calibrations["fourier"]["b"] + if "a" in self.cameraslm.calibrations["fourier"]: + b2 -= np.matmul(M2, self.cameraslm.calibrations["fourier"]["a"]) # Composite transformation (along with xy -> yx). M = cp.array(np.flip(np.flip(np.matmul(M2, M1), axis=0), axis=1)) @@ -240,16 +238,26 @@ def measure(self, basis="ij"): cp.sqrt(self.img_knm, out=self.img_knm) # Target update. - def update_target(self, new_target, reset_weights=False, plot=False): + def update_target(self, new_target_ij, reset_weights=False): + """ + Change the target to something new. This method handles cleaning and normalization. + + Parameters + ---------- + new_target_ij : array_like OR None + New :attr:`target_ij` to optimize towards *in the camera basis*. + should be of the same shape as the camera. + Also updates :attr:`target` using the stored Fourier calibration. + reset_weights : bool + Whether to update the :attr:`weights` to this new :attr:`target`. + """ + self.target_ij = new_target_ij # Transformation order of zero to prevent nan-blurring in MRAF cases. - self.ijcam_to_knmslm(new_target, out=self.target, order=0) + self.ijcam_to_knmslm(new_target_ij, out=self.target, order=0) if reset_weights: self.reset_weights() - if plot: - self.plot_farfield(self.target) - def refine_offset(self, img, basis="kxy"): """ **(NotImplemented)** @@ -275,6 +283,7 @@ def refine_offset(self, img, basis="kxy"): numpy.ndarray Euclidean pixel error in the ``"ij"`` basis for each spot. """ + # Probably local autocorrelation algorithm. raise NotImplementedError() diff --git a/slmsuite/holography/algorithms/_header.py b/slmsuite/holography/algorithms/_header.py index 15ffb74..95cf75f 100644 --- a/slmsuite/holography/algorithms/_header.py +++ b/slmsuite/holography/algorithms/_header.py @@ -27,10 +27,18 @@ cp_gaussian_filter1d = sp_gaussian_filter1d cp_gaussian_filter = sp_gaussian_filter cp_affine_transform = sp_affine_transform - warnings.warn("cupy not installed. Using numpy.") + warnings.warn( + "cupy is not installed; using numpy. Install cupy for faster GPU-based holography." + ) + +try: + import pytorch as torch +except ImportError: + torch = None # Import helper functions from slmsuite.holography import analysis, toolbox +from slmsuite.holography.toolbox import phase as tphase from slmsuite.holography.toolbox.phase import CUDA_KERNELS, _zernike_populate_basis_map from slmsuite.misc.math import REAL_TYPES from slmsuite.misc.files import write_h5, read_h5 @@ -52,8 +60,14 @@ "feedback_exponent": 0.8, }, "WGS-Nogrette": {"feedback": "computational", "feedback_factor": 0.1}, - "WGS-Wu": {"feedback": "computational", "feedback_exponent":1}, - "WGS-tanh": {"feedback": "computational", "feedback_factor": 1, "feedback_exponent":1}, + "WGS-Wu": {"feedback": "computational", "feedback_exponent": 10}, + "WGS-tanh": {"feedback": "computational", "feedback_factor": .5, "feedback_exponent": 10}, + "CG" : { + "feedback": "computational", + "optimizer": "Adam", + "optimizer_kwargs": {"lr": 1e-4}, + "loss": None + } } ALGORITHM_INDEX = {key : i for i, key in enumerate(ALGORITHM_DEFAULTS.keys())} diff --git a/slmsuite/holography/algorithms/_hologram.py b/slmsuite/holography/algorithms/_hologram.py index e69aa3f..699d897 100644 --- a/slmsuite/holography/algorithms/_hologram.py +++ b/slmsuite/holography/algorithms/_hologram.py @@ -102,16 +102,13 @@ class Hologram(_HologramStats): `type promotion `_. Complex datatypes are derived from :attr:`dtype`: - - ``float32`` -> ``complex64`` (default :attr:`dtype`) + - ``float32`` -> ``complex64`` (the default :attr:`dtype`) - ``float64`` -> ``complex128`` ``float16`` is *not* recommended for :attr:`dtype` because ``complex32`` is not implemented by :mod:`numpy`. iter : int Tracks the current iteration number. - method : str - Remembers the name of the last-used optimization method. The method used for each - iteration is stored in ``stats``. flags : dict Helper flags to store custom persistent variables for optimization. These flags are generally changed by passing as a ``kwarg`` to @@ -130,8 +127,9 @@ class Hologram(_HologramStats): - ``"stat_groups"`` : ``list of str`` Stores the values passed to :meth:`~slmsuite.holography.algorithms.Hologram.optimize()`. - - ``"raw_stats"`` : bool - Whether to store raw stats. + - ``"raw_stats"`` : ``bool`` + Whether to store raw stats: the raw image and feedback data for each + iteration. Note that this can be a good amount of data. - ``"blur_ij"`` : ``float`` See :meth:`~slmsuite.holography.algorithms.FeedbackHologram.ijcam_to_knmslm()`. - Other user-defined flags. @@ -154,7 +152,7 @@ class Hologram(_HologramStats): See :meth:`.update_stats()` and :meth:`.plot_stats()`. """ - def __init__(self, target, amp=None, phase=None, slm_shape=None, dtype=np.float32): + def __init__(self, target, amp=None, phase=None, slm_shape=None, dtype=np.float32, **kwargs): r""" Initialize datastructures for optimization. When :mod:`cupy` is enabled, arrays are initialized on the GPU as :mod:`cupy` arrays: @@ -167,8 +165,9 @@ def __init__(self, target, amp=None, phase=None, slm_shape=None, dtype=np.float3 Parameters ---------- target : numpy.ndarray OR cupy.ndarray OR (int, int) - Target to optimize to. The user can also pass a shape in :mod:`numpy` ``(h, - w)`` form, and this constructor will create an empty target of all zeros. + Target to optimize to. + The user can also pass a shape in :mod:`numpy` ``(h, w)`` form, + and this constructor will create an empty target of all zeros. :meth:`.calculate_padded_shape()` can be of particular help for calculating the shape that will produce desired results (in terms of precision, etc). amp : array_like OR None @@ -180,13 +179,15 @@ def __init__(self, target, amp=None, phase=None, slm_shape=None, dtype=np.float3 slm_shape : (int, int) OR slmsuite.hardware.FourierSLM OR slmsuite.hardware.slms.SLM OR None The shape of the near-field of the SLM in :mod:`numpy` `(h, w)` form. Optionally, as a quality of life feature, the user can pass a - :class:`slmsuite.hardware.FourierSLM` or - :class:`slmsuite.hardware.slms.SLM` instead, + :class:`~slmsuite.hardware.FourierSLM` or + :class:`~slmsuite.hardware.slms.SLM` instead, and ``slm_shape`` (and ``amp`` if it is ``None``) are populated from this. If ``None``, tries to use the shape of ``amp`` or ``phase``, but if these are not present, defaults to :attr:`shape` (which is usually determined by ``target``). dtype : type See :attr:`dtype`; type to use for stored arrays. + **kwargs + Passed to :attr:`flags`. """ # 1) Parse inputs # Parse target and create shape. @@ -208,11 +209,13 @@ def __init__(self, target, amp=None, phase=None, slm_shape=None, dtype=np.float3 # 1.5) Determine the shape of the SLM. We have three sources of this shape, which are # optional to pass, but must be self-consistent if passed: - # a) The shape of the nearfield amplitude - # b) The shape of the seed nearfield phase - # c) slm_shape itself (which is set to the shape of a passed SLM, if given). - # If no parameter is passed, these shapes are set to (nan, nan) to prepare for a - # vote (next section). + # + # a) The shape of the nearfield amplitude + # b) The shape of the seed nearfield phase + # c) slm_shape itself (which is set to the shape of a passed SLM, if given). + # + # If no parameter is passed for a given shape, the shape is set to (nan, nan) to + # prepare for a vote (next section). # Option a if amp is None: @@ -295,14 +298,17 @@ def __init__(self, target, amp=None, phase=None, slm_shape=None, dtype=np.float3 self.amp = cp.array(amp, dtype=self.dtype, copy=False) self.amp *= 1 / Hologram._norm(self.amp) - # Initialize near-field phase - self.reset_phase(phase) + # Initialize flags. + self.flags = kwargs # Initialize target. reset() will handle weights. self._update_target(target, reset_weights=False) + # Initialize near-field phase + self.reset_phase(phase) + # Initialize everything else inside reset. - self.reset(reset_phase=False, reset_flags=True) + self.reset(reset_phase=False, reset_flags=False) # Custom GPU kernels for speedy weighting. self._update_weights_generic_cuda_kernel = None @@ -319,9 +325,11 @@ def __init__(self, target, amp=None, phase=None, slm_shape=None, dtype=np.float3 def reset(self, reset_phase=True, reset_flags=False): r""" Resets the hologram to an initial state. Does not restore the preconditioned ``phase`` - that may have been passed to the constructor (as this information is lost upon optimization). + that may have been passed to the constructor (as this information is lost upon + optimization). Instead, phase is randomized if ``reset_phase=True``. Also uses the current ``target`` rather than the ``target`` that may have been - passed to the constructor (e.g. includes current :meth:`.refine_offset()` changes, etc). + passed to the constructor (e.g. includes current + :meth:`~slmsuite.holography.algorithms.SpotHologram.refine_offset()` changes, etc). Parameters ---------- @@ -339,18 +347,54 @@ def reset(self, reset_phase=True, reset_flags=False): # Reset vars. self.iter = 0 - self.method = "" - if reset_flags: - self.flags = {} self.stats = {"method": [], "flags": {}, "stats": {}} + if reset_flags: + self.flags = {"method": ""} # Reset farfield storage. self.amp_ff = None self.phase_ff = None - def reset_phase(self, phase=None): + def _get_quadratic_initial_phase(self, scaling=1): + """ + Analytically guesses a phase pattern (lens, blaze) that will overlap with the target. + """ + target = self.target + if hasattr(target, "get"): + target = self.target.get() + + # Figure out the size of the target in knm space + center_knm = analysis.image_positions(target, nansum=True) # Note this is centered knm space. + + # FUTURE: handle shear. + std_knm = np.sqrt(analysis.image_variances(target, centers=center_knm, nansum=True)[:2, 0]) + std_amp = np.sqrt(analysis.image_variances(self.amp)[:2, 0]) + shape = np.flip(self.shape).astype(float) + + grid = analysis._generate_grid(self.shape[1], self.shape[0], centered=True) + grid[0] /= self.shape[1] + grid[1] /= self.shape[0] + + # Figure out what lens and blaze we should apply to initialize to cover + # the target, based upon the moments we calculated. + return ( + tphase.blaze(grid, center_knm) + + tphase.lens(grid, np.reciprocal(scaling * std_knm * shape / std_amp)) + ) + + def _get_random_phase(self): + if cp == np: # numpy does not support `dtype=` + rng = np.random.default_rng() + return rng.uniform(-np.pi, np.pi, self.slm_shape).astype(self.dtype) + else: + return cp.random.uniform(-np.pi, np.pi, self.slm_shape, dtype=self.dtype) + + def reset_phase(self, phase=None, quadratic_initial_phase=None): r""" - Resets the hologram to a random state or to a provided phase. + Resets the hologram + to a provided phase, + to a random state, + or to a targeted `quadratic phase `_. Parameters ---------- @@ -358,17 +402,40 @@ def reset_phase(self, phase=None): The near-field initial phase. See :attr:`phase`. :attr:`phase` should only be passed if the user wants to precondition the optimization. Of shape :attr:`slm_shape`. + quadratic_initial_phase : bool OR float OR None + We can also precondition the phase analytically (with a lens and blaze) + to roughly the size of the target hologram, according to the first and + second order :meth:`~slmsuite.holography.analysis.image_moments()`. + This quadratic preconditioning is + `thought to help reduce the formation of optical vortices or speckle + `_ + compared to random initialization, as the analytic distribution + is smooth in phase. + If ``None``, looks for ``"quadratic_initial_phase"`` in :attr:`flags`. + If a floating point number is provided, the size of the beam in the + farfield is scaled accordingly. """ + # Parse quadratic_initial_phase + if quadratic_initial_phase is None: + if "quadratic_initial_phase" in self.flags: + quadratic_initial_phase = self.flags["quadratic_initial_phase"] + else: + quadratic_initial_phase = False + # Reset phase to random if no phase is given. if phase is None: - if cp == np: # numpy does not support `dtype=` - rng = np.random.default_rng() - self.phase = rng.uniform(-np.pi, np.pi, self.slm_shape).astype(self.dtype) - else: - self.phase = cp.random.uniform(-np.pi, np.pi, self.slm_shape, dtype=self.dtype) + if quadratic_initial_phase: # Analytic + self.phase = self._get_quadratic_initial_phase(quadratic_initial_phase) + else: # Random + self.phase = self._get_random_phase() else: - # Otherwise, cast as a cp.array with correct type. - self.phase = cp.array(phase, dtype=self.dtype, copy=False) + # Otherwise, cast the passed matrix as a cp.array with correct type. + phase = cp.array(phase, dtype=self.dtype, copy=False) + + if not np.all(np.array(self.slm_shape) == np.array(phase.shape)): + raise ValueError(f"Reset phase of shape {phase.shape} is not of slm_shape {self.slm_shape}") + + self.phase = phase def reset_weights(self): """ @@ -413,7 +480,7 @@ def calculate_padded_shape( ---------- slm_shape : (int, int) OR slmsuite.hardware.FourierSLM The original shape of the SLM in :mod:`numpy` `(h, w)` form. The user can pass a - :class:`slmsuite.hardware.FourierSLM` or :class:`slmsuite.hardware.SLM` instead, + :class:`~slmsuite.hardware.FourierSLM` or :class:`~slmsuite.hardware.SLM` instead, and should pass this when using the ``precision`` parameter. padding_order : int Scales to the ``padding_order`` th larger power of 2. @@ -507,13 +574,17 @@ def _calculate_memory_constrained_shape(device=0, dtype=np.float32): return np.sqrt(num_values_per_array) # User interactions: Changing the target and recovering the nearfield phase and complex farfield. - def _update_target(self, new_target, reset_weights=False, plot=False): + def _update_target(self, new_target, reset_weights=False): """ Change the target to something new. This method handles cleaning and normalization. This method is shelled by :meth:`update_target()` such that it is still accessible in the case that a subclass overwrites :meth:`update_target()`. + Tip + ~~~ + Use :meth:`.plot_farfield()` on :attr:`target`` for visualization. + Parameters ---------- new_target : array_like OR None @@ -521,8 +592,6 @@ def _update_target(self, new_target, reset_weights=False, plot=False): by :class:`SpotHologram`. reset_weights : bool Whether to overwrite ``weights`` with ``target``. - plot : bool - Calls :meth:`.plot_farfield()` on :attr:`target`. """ if new_target is None: self.target = cp.zeros(shape=self.shape, dtype=self.dtype) @@ -540,10 +609,7 @@ def _update_target(self, new_target, reset_weights=False, plot=False): if reset_weights: self.reset_weights() - if plot: - self.plot_farfield(self.target) - - def update_target(self, new_target, reset_weights=False, plot=False): + def update_target(self, new_target, reset_weights=False): """ Change the target to something new. This method handles cleaning and normalization. @@ -555,10 +621,8 @@ def update_target(self, new_target, reset_weights=False, plot=False): be used by a user). reset_weights : bool Whether to update the :attr:`weights` to this new :attr:`target`. - plot : bool - Calls :meth:`.plot_farfield()` on :attr:`target`. """ - self._update_target(new_target=new_target, reset_weights=reset_weights, plot=plot) + self._update_target(new_target=new_target, reset_weights=reset_weights) def extract_phase(self): r""" @@ -583,7 +647,6 @@ def extract_farfield(self, affine=None, get=True): ---------- affine : dict Affine transformation to apply to far-field data (in the form of - :attr:`~slmsuite.hardware.cameraslms.FourierSLM.fourier_calibration`). get : bool Whether or not to convert the cupy array to a numpy array if cupy is used. This is ignored if numpy is used. @@ -625,6 +688,46 @@ def extract_farfield(self, affine=None, get=True): ) return farfield + def _populate_results(self, nearfield=None): + """ + Helper function to populate: + - farfield + - amp_ff + - phase_ff + + From the data in: + - amp + - phase + + Pass a nearfield complex matrix of self.shape shape to avoid memory reallocation. + """ + (i0, i1, i2, i3) = toolbox.unpad(self.shape, self.slm_shape) + + if nearfield is None: + nearfield = cp.zeros(self.shape, dtype=self.dtype_complex) + else: + nearfield.fill(0) + + nearfield[i0:i1, i2:i3] = self.amp * cp.exp(1j * self.phase) + + self.farfield = self._nearfield2farfield(nearfield, farfield_out=self.farfield) + self.amp_ff = cp.abs(self.farfield) + self.phase_ff = cp.angle(self.farfield) + + def _nearfield2farfield(self, nearfield, farfield_out=None): + """ + Maps the nearfield to the farfield by a discrete Fourier transform. + This function is overloaded by subclasses. + """ + return cp.fft.fftshift(cp.fft.fft2(cp.fft.fftshift(nearfield), norm="ortho")) + + def _farfield2nearfield(self, farfield, nearfield_out=None): + """ + Maps the farfield to the nearfield by a discrete Fourier transform. + This function is overloaded by subclasses. + """ + return cp.fft.ifftshift(cp.fft.ifft2(cp.fft.ifftshift(farfield), norm="ortho")) + # Core optimization function. def optimize( self, @@ -672,7 +775,7 @@ def optimize( feedback (measured) amplitudes, and :math:`p` is the power passed as ``"feedback_exponent"`` in :attr:`~slmsuite.holography.algorithms.Hologram.flags` (see ``kwargs``). - The power :math:`p` defaults to .9 if not passed. In general, smaller + The power :math:`p` defaults to .8 if not passed. In general, smaller :math:`p` will lead to slower yet more stable optimization. - ``'WGS-Kim'`` @@ -700,9 +803,32 @@ def optimize( - ``'WGS-Wu'`` - `exponential `_ + TODO `exponential `_ + + .. math:: \mathcal{W} = \mathcal{W}\exp\left( p (\mathcal{T} - \mathcal{F}) \right) + + The speed of correction is controlled by :math:`p`, + the power passed as ``"feedback_exponent"``. + + - ``'WGS-tanh'`` + + TODO Hyperbolic tangent + + .. math:: \mathcal{W} = \mathcal{W}\left[1 + f\text{tanh}\left( p (\mathcal{T} - \mathcal{F}) \right) \right] + + This weighting limits each update to a relative change of :math:`\pm f`, + which is passed as ``"feedback_factor"``. + The speed of correction is controlled by :math:`p`, + the power passed as ``"feedback_exponent"``. - .. math:: \mathcal{W} = \mathcal{W}\exp\left( f (\mathcal{T} - \mathcal{F}) \right) + - Conjugate Gradient (CG) phase retrieval. **(Not Finished)** + + - ``'CG'`` + + TODO + + Uses :mod:`pytorch`-:mod:`cupy` + `interoperability `_. - The option for `Mixed Region Amplitude Freedom (MRAF) `_ feedback. In standard @@ -808,7 +934,7 @@ def optimize( "unrecognized method '{}'.\n" "Valid methods include {}".format(method, methods) ) - self.method = method + self.flags["method"] = method # 1) Parse flags: # 1.1) Set defaults if not already set. @@ -842,9 +968,7 @@ def optimize( # 1.4) Print the flags if verbose. if verbose > 1: print( - "Optimizing with '{}' using the following method-specific flags:".format( - self.method - ) + f"Optimizing with '{method}' using the following method-specific flags:" ) pprint.pprint( { @@ -865,8 +989,12 @@ def optimize( # 3) Switch between optimization methods (currently only GS- or WGS-type is supported). if "GS" in method: self.optimize_gs(iterations, callback) + elif "CG" in method: + self.optimize_cg(iterations, callback) + else: + raise ValueError(f"Unsupported optimization method '{method}'") - # Optimization methods (currently only GS- or WGS-type is supported). + # GS- or WGS-type optimization. def optimize_gs(self, iterations, callback): """ GPU-accelerated Gerchberg-Saxton (GS) iterative phase retrieval. @@ -901,9 +1029,10 @@ def optimize_gs(self, iterations, callback): callback : callable OR None See :meth:`.optimize()`. """ - # Proxy to initialize nearfield with the correct shape and correct (complex) type. + # Initialize nearfield and farfield. nearfield = cp.zeros(self.shape, dtype=self.dtype_complex) - self.farfield = cp.zeros(self.target.shape, dtype=self.dtype_complex) # Use target.shape for FreeSpotHologram cases. + # Use self.target.shape instead of self.shape to account for CompressedSpotHologram cases. + self.farfield = cp.zeros(self.target.shape, dtype=self.dtype_complex) # Precompute MRAF helper variables. mraf_variables = self._mraf_helper_routines() @@ -931,7 +1060,7 @@ def optimize_gs(self, iterations, callback): if callback(self): break - # Evaluate method-specific routines, stats, etc. + # Evaluate method-specific routines, stats, etc. This includes camera feedback/etc. # If you want to add new functionality to GS, do so here to keep the main loop clean. self._gs_farfield_routines(self.farfield, mraf_variables) @@ -949,26 +1078,8 @@ def optimize_gs(self, iterations, callback): # Increment iteration. self.iter += 1 - # Update the final far-field - nearfield.fill(0) - nearfield[i0:i1, i2:i3] = self.amp * cp.exp(1j * self.phase) - self.farfield = self._nearfield2farfield(nearfield, farfield_out=self.farfield) - self.amp_ff = cp.abs(self.farfield) - self.phase_ff = cp.angle(self.farfield) - - def _nearfield2farfield(self, nearfield, farfield_out=None): - """ - Maps the nearfield to the farfield by a discrete Fourier transform. - This function is overloaded by subclasses. - """ - return cp.fft.fftshift(cp.fft.fft2(cp.fft.fftshift(nearfield), norm="ortho")) - - def _farfield2nearfield(self, farfield, nearfield_out=None): - """ - Maps the farfield to the nearfield by a discrete Fourier transform. - This function is overloaded by subclasses. - """ - return cp.fft.ifftshift(cp.fft.ifft2(cp.fft.ifftshift(farfield), norm="ortho")) + # Update the final far-field using phase and amp. + self._populate_results() def _mraf_helper_routines(self): # MRAF helper variables @@ -1035,11 +1146,11 @@ def _gs_farfield_routines(self, farfield, mraf_variables): self.update_stats(self.flags["stat_groups"]) # Weight, if desired. - if "WGS" in self.method: + if "WGS" in self.flags["method"]: self._update_weights() # Decide whether to fix phase. - if "Kim" in self.method: + if "Kim" in self.flags["method"]: was_not_fixed = not self.flags["fixed_phase"] # Enable based on efficiency. @@ -1047,7 +1158,8 @@ def _gs_farfield_routines(self, farfield, mraf_variables): stats = self.stats["stats"] groups = tuple(stats.keys()) - assert len(stats) > 0, "Must track statistics to fix phase based on efficiency!" + if len(stats) == 0: + raise ValueError("Must track statistics to fix phase based on efficiency!") eff = stats[groups[-1]]["efficiency"][self.iter] if eff > self.flags["fix_phase_efficiency"]: @@ -1083,7 +1195,7 @@ def _gs_farfield_routines(self, farfield, mraf_variables): # cp.multiply(farfield, self.weights, out=farfield) # cp.nan_to_num(farfield, copy=False, nan=0) - if not ("fixed_phase" in self.flags and self.flags["fixed_phase"]): + if not ("fixed_phase" in self.flags and self.flags["fixed_phase"]) or self.phase_ff is None: self.phase_ff = cp.arctan2(farfield.imag, farfield.real, out=self.phase_ff) cp.exp(1j * self.phase_ff, out=farfield) @@ -1134,9 +1246,113 @@ def _gs_farfield_routines(self, farfield, mraf_variables): cp.multiply(farfield, self.weights, _where=signal_region, out=farfield) if mraf_factor is not None: cp.multiply(farfield, mraf_factor, _where=noise_region, out=farfield) + # Conjugate gradient optimization. + def optimize_cg(self, iterations, callback): + """ + **(Not Finished)** + + Conjugate Gradient (CG) iterative phase retrieval. + + Solves the "phase problem": approximates the near-field phase that + transforms a known near-field source amplitude to a desired far-field + target amplitude. + + Caution + ~~~~~~~ + This function should be called through :meth:`.optimize()` and not called + directly. It is left as a public function exposed in documentation to clarify + how the internals of :meth:`.optimize()` work. + + Parameters + ---------- + iterations : iterable + Number of loop iterations to run. Is an iterable to pass a :mod:`tqdm` iterable. + callback : callable OR None + See :meth:`.optimize()`. + """ + if torch is None: + raise ValueError("pytorch is required for conjugate gradient optimization.") + + # Initialize nearfield and farfield. + nearfield = cp.zeros(self.shape, dtype=self.dtype_complex) + # Use self.target.shape instead of self.shape to account for CompressedSpotHologram cases. + self.farfield = cp.zeros(self.target.shape, dtype=self.dtype_complex) + + # Convert variables to torch with **zero-copy** cupy interoperability. + # We need torch to handle gradient calculation. + amp_torch = Hologram._get_torch_tensor_from_cupy(self.amp) + target_torch = Hologram._get_torch_tensor_from_cupy(self.target) + + phase_torch = Hologram._get_torch_tensor_from_cupy(self.phase) + nearfield_torch = Hologram._get_torch_tensor_from_cupy(nearfield) + farfield_torch = Hologram._get_torch_tensor_from_cupy(self.farfield) + + # We're going to use phase for optimization, so we need to record gradients. + phase_torch.requires_grad_(True) + nearfield_torch.requires_grad_(True) + farfield_torch.requires_grad_(True) + + # Parse arguments. + loss = self.flags["loss"] + if loss is None: + loss = torch.nn.MSELoss() + + # Create the optimizer. + try: + optim_class = getattr(torch.optim, self.flags["optimizer"]) + except: + raise ValueError(f"{self.flags["optimizer"]} is not a valid torch optimizer") + + self.optimizer = optim_class([phase_torch], **self.flags["optimizer_kwargs"]) + + # Helper variables for speeding up source phase and amplitude fixing. + (i0, i1, i2, i3) = toolbox.unpad(self.shape, self.slm_shape) + + for _ in iterations: + # (A) Nearfield -> farfield + # Fix the relevant part of the nearfield amplitude to the source amplitude. + # Everything else is zero because power outside the SLM is assumed unreflected. + # This is optimized for when shape is much larger than slm_shape. + nearfield_torch.fill_(0) + self.optimizer.zero_grad() + nearfield_torch[i0:i1, i2:i3] = amp_torch * torch.exp(1j * phase_torch) + + # FFT to move to the farfield. TODO: Torchify; current implementation will + # not track gradients. + self.farfield = self._nearfield2farfield(nearfield, farfield_out=self.farfield) + + # (B) Midloop caching and prep + # Before callback(), cleanup such that it can access updated amp_ff and images. + self._midloop_cleaning(self.farfield) + + # Run step function if present and check termination conditions. + if callback is not None: + if callback(self): + break + + # (C) Evaluate loss + result = loss(farfield_torch, target_torch) + result.backward(retain_graph=True) + + # (D) Finally, step the optimization according to the gradients calculated. + self.optimizer.step() + + # Increment iteration. + self.iter += 1 + + # Update the final far-field using phase and amp. + self._populate_results() + + @staticmethod + def _get_torch_tensor_from_cupy(array): + if cp == np: + return torch.as_tensor(array, device='cuda') + else: + return torch.as_tensor(array, device='cpu') + # Weighting functions. def _update_weights_generic( - self, weight_amp, feedback_amp, target_amp=None, xp=cp, nan_checks=True + self, weight_amp, feedback_amp, target_amp, xp=cp, nan_checks=True ): """ Helper function to process weight feedback according to the chosen weighting method. @@ -1148,15 +1364,13 @@ def _update_weights_generic( Parameters ---------- weight_amp : numpy.ndarray OR cupy.ndarray - A :class:`~slmsuite.holography.SpotArray` instance containing locations - where the feedback weight should be calculated. + Weights to update. feedback_amp : numpy.ndarray OR cupy.ndarray Measured or result amplitudes corresponding to ``weight_amp``. Should be the same size as ``weight_amp``. target_amp : numpy.ndarray OR cupy.ndarray OR None Necessary in the case where ``target_amp`` is not uniform, such that the weighting can - properly be applied to bring the feedback closer to the target. If ``None``, is assumed - to be uniform. Should be the same size as ``weight_amp``. + properly be applied to bring the feedback closer to the target. xp : module This function is used by both :mod:`cupy` and :mod:`numpy`, so we have the option for either. Defaults to :mod:`cupy`. @@ -1174,18 +1388,20 @@ def _update_weights_generic( return self._update_weights_generic_cuda(weight_amp, feedback_amp, target_amp) def _update_weights_generic_cupy( - self, weight_amp, feedback_amp, target_amp=None, xp=cp, nan_checks=True + self, weight_amp, feedback_amp, target_amp, xp=cp, nan_checks=True ): - assert self.method[:4] == "WGS-", "For now, assume weighting is for WGS." - method = self.method[4:].lower() - - # Parse feedback_amp - if target_amp is None: # Uniform - feedback_corrected = xp.array(feedback_amp, copy=True, dtype=self.dtype) - else: # Non-uniform - feedback_corrected = xp.array(feedback_amp, copy=True, dtype=self.dtype) - feedback_corrected *= 1 / Hologram._norm(feedback_corrected, xp=xp) - + method = self.flags["method"].lower() + if method[:4] != "wgs-": + raise ValueError("Weighting is only for WGS.") + method = method[4:] + + feedback_corrected = xp.array(feedback_amp, copy=True, dtype=self.dtype) + feedback_corrected *= 1 / Hologram._norm(feedback_corrected, xp=xp) + + if ("wu" in method or "tanh" in method): # Additive + feedback_corrected *= -self.flags["feedback_exponent"] + feedback_corrected += xp.array(target_amp, copy=False) + else: # Multiplicative xp.divide(feedback_corrected, xp.array(target_amp, copy=False), out=feedback_corrected) if nan_checks: @@ -1205,13 +1421,14 @@ def _update_weights_generic_cupy( feedback_corrected *= -self.flags["feedback_factor"] feedback_corrected += 1 xp.reciprocal(feedback_corrected, out=feedback_corrected) - # elif "wu" in method: - # feedback = np.exp(self.flags["feedback_exponent"] * (target - feedback)) - # elif "tanh" in method: - # feedback = self.flags["feedback_factor"] * np.tanh(self.flags["feedback_exponent"] * (target - feedback)) + elif "wu" in method: + feedback_corrected = np.exp(feedback_corrected) + elif "tanh" in method: + feedback_corrected = self.flags["feedback_factor"] * np.tanh(feedback_corrected) + feedback_corrected += 1 else: - raise RuntimeError( - "Method " "{}" " not recognized by Hologram.optimize()".format(self.method) + raise ValueError( + f"Method '{self.flags["method"]}' not recognized by Hologram.optimize()" ) if nan_checks: @@ -1229,22 +1446,18 @@ def _update_weights_generic_cupy( return weight_amp - def _update_weights_generic_cuda(self, weight_amp, feedback_amp, target_amp=None): - method = ALGORITHM_INDEX[self.method] + def _update_weights_generic_cuda(self, weight_amp, feedback_amp, target_amp): + N = weight_amp.size - if target_amp is None: # Uniform - feedback_norm = 0 - target_amp = 0 - else: - feedback_norm = Hologram._norm(feedback_amp, xp=cp) + feedback_amp = cp.array(feedback_amp, copy=False) + feedback_norm = Hologram._norm(feedback_amp, xp=cp) - N = weight_amp.size - method = ALGORITHM_INDEX[self.method] + method = ALGORITHM_INDEX[self.flags["method"]] threads_per_block = int( self._update_weights_generic_cuda_kernel.max_threads_per_block ) - blocks = N // threads_per_block + blocks = N // threads_per_block + 1 # Call the RawKernel. self._update_weights_generic_cuda_kernel( @@ -1257,8 +1470,8 @@ def _update_weights_generic_cuda(self, weight_amp, feedback_amp, target_amp=None N, method, feedback_norm, - self.flags.pop("feedback_exponent", 0), # TODO: fix - self.flags.pop("feedback_factor", 0) + self.flags.pop("feedback_exponent", 1), # TODO: fix + self.flags.pop("feedback_factor", 1) ) ) diff --git a/slmsuite/holography/algorithms/_spots.py b/slmsuite/holography/algorithms/_spots.py index ba48fbe..8a58e6c 100644 --- a/slmsuite/holography/algorithms/_spots.py +++ b/slmsuite/holography/algorithms/_spots.py @@ -242,9 +242,9 @@ def __init__( The index ``-1`` (outside Zernike indexing) is used as a special case to add a vortex waveplate with amplitude :math:`2\pi` to the system - (see :meth:`slmsuite.holography.toolbox.phase.laguerre_gaussian()`). + (see :meth:`~slmsuite.holography.toolbox.phase.laguerre_gaussian()`). zernike_center_shift : array_like OR None - + TODO **kwargs Passed to :meth:`.FeedbackHologram.__init__()`. """ @@ -774,18 +774,18 @@ def expand_kernel(kernel, farfield, out): return nearfield_out # Target update. - def update_target(self, new_target, reset_weights=False, plot=False): + def update_target(self, new_target, reset_weights=False): """ Change the target to something new. This method handles cleaning and normalization. Parameters ---------- new_target : array_like OR None + A list with ``N`` elements corresponding to the target intensities of each + of the ``N`` spots. If ``None``, sets the target spot amplitudes to be uniform. reset_weights : bool Whether to overwrite ``weights`` with ``target``. - plot : bool - Calls :meth:`.plot_farfield()` on :attr:`target`. """ if new_target is None: # Default to even power on all spots. @@ -806,9 +806,6 @@ def update_target(self, new_target, reset_weights=False, plot=False): if reset_weights: self.reset_weights() - if plot: - raise NotImplementedError() - # Weighting and stats. def _update_weights(self): """ @@ -1075,7 +1072,7 @@ def __init__( self.spot_knm, "knm", "kxy", cameraslm.slm, shape ) - if cameraslm.fourier_calibration is not None: + if "fourier" in cameraslm.calibrations: self.spot_ij = cameraslm.kxyslm_to_ijcam(self.spot_kxy) else: self.spot_ij = None @@ -1092,8 +1089,8 @@ def __init__( self.spot_kxy = vectors - if hasattr(cameraslm, "fourier_calibration"): - if cameraslm.fourier_calibration is not None: + if hasattr(cameraslm, "calibrations"): + if "fourier" in cameraslm.calibrations: self.spot_ij = cameraslm.kxyslm_to_ijcam(vectors) # This is okay for non-feedback GS, so we don't error. else: @@ -1104,7 +1101,7 @@ def __init__( ) elif basis == "ij": # Pixel on the camera. assert cameraslm is not None, "We need an cameraslm to interpret ij." - assert cameraslm.fourier_calibration is not None, ( + assert "fourier" in cameraslm.calibrations, ( "We need an cameraslm with " "fourier-calibrated kxyslm_to_ijcam and ijcam_to_kxyslm transforms " "to interpret ij." @@ -1211,7 +1208,7 @@ def __init__( self.null_radius_knm = int(np.ceil(self.null_radius_knm)) - # Initialize target/etc. + # Initialize target/etc. Note that this passes through FeedbackHologram. super().__init__(shape, target_ij=None, cameraslm=cameraslm, **kwargs) # Parse null_region after __init__ @@ -1320,7 +1317,7 @@ def make_rectangular_array( assert "cameraslm" in kwargs, "We need an cameraslm to interpret ij." cameraslm = kwargs["cameraslm"] assert cameraslm is not None, "We need an cameraslm to interpret ij." - assert cameraslm.fourier_calibration is not None, ( + assert "fourier" in cameraslm.calibrations, ( "We need an cameraslm with " "fourier-calibrated kxyslm_to_ijcam and ijcam_to_kxyslm transforms " "to interpret ij." @@ -1348,7 +1345,7 @@ def make_rectangular_array( # Return a new SpotHologram. return SpotHologram(shape, vectors, basis=basis, spot_amp=None, **kwargs) - def _update_target_spots(self, reset_weights=False, plot=False): + def _update_target_spots(self, reset_weights=False): """ Wrapped by :meth:`SpotHologram.update_target()`. """ @@ -1365,7 +1362,7 @@ def _update_target_spots(self, reset_weights=False, plot=False): self.shape, ) - if self.cameraslm.fourier_calibration is not None: + if "fourier" in self.cameraslm.calibrations: self.spot_ij_rounded = self.cameraslm.kxyslm_to_ijcam(self.spot_kxy_rounded) else: self.spot_ij_rounded = None @@ -1406,17 +1403,14 @@ def _update_target_spots(self, reset_weights=False, plot=False): if reset_weights: self.reset_weights() - if plot: - self.plot_farfield(self.target) - def update_target(self, reset_weights=False, plot=False): """ From the spot locations stored in :attr:`spot_knm`, update the target pattern. Note ~~~~ - If there's a cameraslm, updates the :attr:`spot_ij_rounded` attribute - corresponding to where pixels in the k-space where actually placed (due to rounding + If there's a ``cameraslm``, updates the :attr:`spot_ij_rounded` attribute + corresponding to where pixels in the :math:`k`-space where actually placed (due to rounding to integers, stored in :attr:`spot_knm_rounded`), rather the idealized floats :attr:`spot_knm`. @@ -1424,17 +1418,14 @@ def update_target(self, reset_weights=False, plot=False): ~~~~ The :attr:`target` and :attr:`weights` matrices are modified in-place for speed, unlike :class:`.Hologram` or :class:`.FeedbackHologram` which make new matrices. - This is because spot positions are expected to be corrected using :meth:`correct_spots()`. + This is because spot positions are expected to be corrected using :meth:`refine_offsets()`. Parameters ---------- reset_weights : bool - Whether to rest the :attr:`weights` to this new :attr:`target`. - plot : bool - Whether to enable debug plotting to see the positions of the spots relative to the - shape of the camera and slm. + Whether to reset the :attr:`weights` to this new :attr:`target`. """ - self._update_target_spots(reset_weights=reset_weights, plot=plot) + self._update_target_spots(reset_weights=reset_weights) # Weighting and stats. def _update_weights(self): diff --git a/slmsuite/holography/algorithms/_stats.py b/slmsuite/holography/algorithms/_stats.py index b7920da..f477a15 100644 --- a/slmsuite/holography/algorithms/_stats.py +++ b/slmsuite/holography/algorithms/_stats.py @@ -23,7 +23,7 @@ def _calculate_stats( Target of holography. xp : module This function is used by both :mod:`cupy` and :mod:`numpy`, so we have the option - for either. Defaults to :mod:`cupy`. + to use either. Defaults to :mod:`cupy`. efficiency_compensation : bool Whether to scale the ``feedback`` based on the overlap with the ``target``. This is more accurate for images, but less accurate for SpotHolograms. @@ -35,7 +35,7 @@ def _calculate_stats( raw : bool Passes the ``"raw_stats"`` flag. If ``True``, stores the raw feedback and raw feedback-target ratio for each pixel or spot instead of - only the derived statistics. + only the derived and compressed statistics. """ # Downgrade to numpy if necessary if xp == np and (hasattr(feedback_amp, "get") or hasattr(target_amp, "get")): @@ -79,7 +79,6 @@ def _calculate_stats( # Make some helper lists; ignoring power where target is zero. mask = xp.logical_and(target_pwr != 0, xp.logical_not(xp.isnan(target_pwr))) - # mask = xp.nonzero(target_pwr) feedback_pwr_masked = feedback_pwr[mask] target_pwr_masked = target_pwr[mask] @@ -143,7 +142,7 @@ def _update_stats_dictionary(self, stats): if diff > 0: # Extend methods self.stats["method"].extend(["" for _ in range(diff)]) M = self.iter + 1 - self.stats["method"][self.iter] = self.method # Update method + self.stats["method"][self.iter] = self.flags["method"] # Update method # Update flags flaglist = set(self.flags.keys()).union(set(self.stats["flags"].keys())) @@ -240,26 +239,31 @@ def export_stats(self, file_path, include_state=True): interfaces), but rather to provide extra information for debugging. """ # Save attributes, converting to numpy when necessary. + to_save = {} + if include_state: - to_save = { - "slm_shape": self.slm_shape, - "phase": self.phase, - "amp": self.amp, - "shape": self.shape, - "target": self.target, - "weights": self.weights, - "phase_ff": self.phase_ff, - "iter": self.iter, - "method": self.method, - "flags": self.flags, - } - - for key in to_save.keys(): - if hasattr(to_save[key], "get") and not isinstance(to_save[key], dict): - to_save[key] = to_save[key].get() - else: + to_save_keys = [ + "slm_shape", + "phase", + "amp", + "shape", + "target", + "weights", + "phase_ff", + "iter", + "method", + "flags", + ] to_save = {} + for key in to_save_keys: + value = getattr(self, key) + + if hasattr(value, "get") and not isinstance(value, dict): + to_save[key] = value.get() + else: + to_save[key] = value + # Save stats. to_save["stats"] = self.stats @@ -487,7 +491,7 @@ def plot_farfield( limits = self._compute_limits(self.target.get(), limit_padding=limit_padding) if len(title) == 0: - title = "FF Amp" + title = "Farfield Amplitude" # Interpret source and convert to numpy for plotting. isphase = "phase" in title.lower() @@ -599,15 +603,15 @@ def rebase(ax, img, to_units): if i == 0: ax.set_ylabel(toolbox.BLAZE_LABELS[units][1]) - # If cam_points is defined (i.e. is a FeedbackHologram or subclass), + # If _cam_points is defined (i.e. is a FeedbackHologram or subclass), # plot a yellow rectangle for the extents of the camera - if hasattr(self, "cam_points") and self.cam_points is not None: + if hasattr(self, "_cam_points") and self._cam_points is not None: # Check to see if the camera extends outside of knm space. plot_slm_fov = ( - np.any(self.cam_points[0, :4] < 0) - or np.any(self.cam_points[1, :4] < 0) - or np.any(self.cam_points[0, :4] >= npsource.shape[1]) - or np.any(self.cam_points[1, :4] >= npsource.shape[1]) + np.any(self._cam_points[0, :4] < 0) + or np.any(self._cam_points[1, :4] < 0) + or np.any(self._cam_points[0, :4] >= npsource.shape[1]) + or np.any(self._cam_points[1, :4] >= npsource.shape[1]) ) # If so, plot a labeled green rectangle to show the extents of knm space. @@ -631,23 +635,23 @@ def rebase(ax, img, to_units): va="top", ) - # Convert cam_points to knm. + # Convert _cam_points to knm. if units == "knm": - cam_points = self.cam_points + _cam_points = self._cam_points else: - cam_points = toolbox.convert_vector( - self.cam_points, from_units="knm", to_units=units, slm=slm, shape=npsource.shape + _cam_points = toolbox.convert_vector( + self._cam_points, from_units="knm", to_units=units, slm=slm, shape=npsource.shape ) # Plot the labeled yellow rectangle representing the camera. axs[0].plot( - cam_points[0], - cam_points[1], + _cam_points[0], + _cam_points[1], c="y", ) axs[0].annotate( "Camera FoV", - (np.mean(cam_points[0, :4]), np.max(cam_points[1, :4])), + (np.mean(_cam_points[0, :4]), np.max(_cam_points[1, :4])), c="y", size="small", ha="center", @@ -656,22 +660,22 @@ def rebase(ax, img, to_units): # Determine sensible limits of the field of view. if plot_slm_fov: - dx = np.max(cam_points[0]) - np.min(cam_points[0]) - dy = np.max(cam_points[1]) - np.min(cam_points[1]) + dx = np.max(_cam_points[0]) - np.min(_cam_points[0]) + dy = np.max(_cam_points[1]) - np.min(_cam_points[1]) else: dx = dy = 0 ext = full.get_extent() axs[0].set_xlim( [ - min(ext[0], np.min(cam_points[0]) - dx / 10), - max(ext[1], np.max(cam_points[0]) + dx / 10), + min(ext[0], np.min(_cam_points[0]) - dx / 10), + max(ext[1], np.max(_cam_points[0]) + dx / 10), ] ) axs[0].set_ylim( [ - max(ext[2], np.max(cam_points[1]) + dy / 10), - min(ext[3], np.min(cam_points[1]) - dy / 10), + max(ext[2], np.max(_cam_points[1]) + dy / 10), + min(ext[3], np.min(_cam_points[1]) - dy / 10), ] ) diff --git a/slmsuite/holography/analysis/__init__.py b/slmsuite/holography/analysis/__init__.py index c59528c..6622b2a 100644 --- a/slmsuite/holography/analysis/__init__.py +++ b/slmsuite/holography/analysis/__init__.py @@ -85,7 +85,9 @@ def take( Defaults to ``True``. integrate : bool If ``True``, the spatial dimension are integrated (summed), yielding a result of the - same length as the number of vectors. Defaults to ``False``. + same length as the number of vectors. Forces floating point datatype before the + summation is done, as integer data (especially for cameras near saturation) can overflow. + Defaults to ``False``. clip : bool Whether to allow out-of-range integration regions. ``True`` allows regions outside the valid area, setting the invalid region to ``np.nan`` @@ -184,7 +186,7 @@ def take( pass if integrate: # Sum over the integration axis. - return xp.squeeze(xp.sum(result, axis=-1)) + return xp.squeeze(xp.sum(result.astype(float), axis=-1)) else: # Reshape the integration axis. return xp.reshape(result, (vectors.shape[1], size[1], size[0])) @@ -238,6 +240,11 @@ def image_remove_field(images, deviations=1, out=None): computed uniquely for each image, or the median of each image if ``deviations`` is ``None``. This is equivalent to background subtraction. + Important + ~~~~~~~~~ + If a stack of images is provided, field removal is done individually on each image. + Field removal is not done in aggregate. + Parameters ---------- images : numpy.ndarray @@ -299,7 +306,7 @@ def image_moment(images, moment=(1, 0), centers=(0, 0), grid=None, normalize=Tru r""" Computes the given `moment `_ :math:`M_{m_xm_y}` for a stack of images. - This involves integrating each image against polynomial trial functions: + This involves discretely integrating each image against polynomial trial functions: .. math:: M_{m_xm_y} = \frac{ \int_{-w_x/2}^{+w_x/2} dx \, (x-c_x)^{m_x} \int_{-w_y/2}^{+w_y/2} dy \, (y-c_y)^{m_y} @@ -352,9 +359,13 @@ def image_moment(images, moment=(1, 0), centers=(0, 0), grid=None, normalize=Tru ``float`` or an anisotropic ``(float, float)``. This corresponds to the pixel's :math:`\Delta x`, :math:`\Delta y`. - Providing lists of length ``w`` and ``h`` as a tuple as the grid dimension. - - Providing full grids of shape ``(w, h)`` in each direction. Note that this - case is the most general, and can lead to a rotated grid if a transformed - grid is provided. + + Tip + ~~~ + If the user pre-allocates and reuses these lists, this case has best performance. + - Providing two full grids of shape ``(h, w)``, one for each direction. + Note that this case is the most general, and can lead to a rotated grid if a + transformed grid is provided. normalize : bool Whether to normalize ``images``. @@ -434,12 +445,13 @@ def image_moment(images, moment=(1, 0), centers=(0, 0), grid=None, normalize=Tru else: x_grid, y_grid = grid - if len(np.shape(x_grid)) == 2: - # 2D grids. + x_grid = np.squeeze(x_grid) + y_grid = np.squeeze(y_grid) + + if len(np.shape(x_grid)) == 2: # 2D grids. x_grid = np.reshape(x_grid, (1, w_y, w_x)) - c_x y_grid = np.reshape(y_grid, (1, w_y, w_x)) - c_y - elif len(np.shape(x_grid)) == 1: - # 1D grids. + elif len(np.shape(x_grid)) == 1: # 1D grids. x_grid = np.reshape(x_grid, (1, 1, w_x)) - c_x y_grid = np.reshape(y_grid, (1, w_y, 1)) - c_y elif len(np.shape(x_grid)) == 3: @@ -448,21 +460,15 @@ def image_moment(images, moment=(1, 0), centers=(0, 0), grid=None, normalize=Tru raise ValueError(f"Could not parse grid of shape {x_grid.shape}") # Don't modify original memory. - if moment[0] > 0: x_grid = x_grid.copy() - if moment[0] > 1: x_grid = np.power(x_grid, moment[0], out=x_grid) - - if moment[1] > 0: y_grid = y_grid.copy() - if moment[1] > 1: y_grid = np.power(y_grid, moment[1], out=y_grid) + if moment[0] > 1: x_grid = np.power(x_grid, moment[0]) + if moment[1] > 1: y_grid = np.power(y_grid, moment[1]) if moment[1] == 0: # Only-x case. - x_grid *= reciprocal - return np_sum(images * x_grid, axis=(1, 2), keepdims=False) + return np_sum(images * x_grid * reciprocal, axis=(1, 2), keepdims=False) elif moment[0] == 0: # Only-y case. - y_grid *= reciprocal - return np_sum(images * y_grid, axis=(1, 2), keepdims=False) + return np_sum(images * y_grid * reciprocal, axis=(1, 2), keepdims=False) else: # Shear case. - x_grid *= reciprocal - return np_sum(images * x_grid * y_grid, axis=(1, 2), keepdims=False) + return np_sum(images * x_grid * y_grid * reciprocal, axis=(1, 2), keepdims=False) def image_normalization(images, nansum=False): @@ -717,6 +723,31 @@ def image_ellipticity(variances): return 1 - (eig_minus / eig_plus) +def image_areas(variances): + r""" + Given the output of :meth:`image_variances()`, + return a measure of spot area for each moment triplet. + The output of :meth:`image_variances()` contains the moments :math:`M_{20}`, + :math:`M_{02}`, and :math:`M_{11}`. We return the determinant :math:`|M|` which is a + proxy for spot area. + + Parameters + ---------- + variances : numpy.ndarray + The output of :meth:`image_variances()`. Shape ``(3, image_count)``. + + Returns + ------- + numpy.ndarray + Array of areas for the given moments in an array of shape ``(image_count,)``. + """ + m20 = variances[0, :] + m02 = variances[1, :] + m11 = variances[2, :] + + return m20 * m02 - m11 * m11 + + def image_ellipticity_angle(variances): r""" Given the output of :meth:`image_variances()`, @@ -756,6 +787,10 @@ def image_ellipticity_angle(variances): return np.arctan2(eig_plus - m02, m11, where=m11 != 0, out=np.zeros_like(m11)) +def batch_fit(y, x, function, guess, plot=False): + pass + + def image_fit(images, grid=None, function=gaussian2d, guess=None, plot=False): """ Fit each image in a stack of images to a 2D ``function``. @@ -771,14 +806,16 @@ def image_fit(images, grid=None, function=gaussian2d, guess=None, plot=False): function : lambda ((float, float), ... ) -> float Some fitfunction which accepts ``(x,y)`` coordinates as first argument. Defaults to :meth:`~slmsuite.misc.fitfunctions.gaussian2d()`. - guess : None OR numpy.ndarray (``image_count``, ``parameter_count``) - - If ``guess`` is ``None``, will construct a guess based on the ``function`` passed. + guess : None OR True OR numpy.ndarray (``image_count``, ``parameter_count``) + - If ``guess`` is ``None`` or ``True``, will construct a guess based on the ``function`` passed. Functions for which guesses are implemented include: - :meth:`~slmsuite.misc.fitfunctions.gaussian2d()` - If ``guess`` is ``None`` and ``function`` does not have a guess - implemented, no guess will be provided to the optimizer. + implemented, no guess will be provided to the optimizer and the user will be warned. + - If ``guess`` is ``True`` and ``function`` does not have a guess + implemented, an error will be raised. - If ``guess`` is a ``numpy.ndarray``, a column of the array will be provided to the optimizer as a guess for the fit parameters for each image. plot : bool @@ -793,10 +830,10 @@ def image_fit(images, grid=None, function=gaussian2d, guess=None, plot=False): A matrix with the fit results. The first row contains the rsquared quality of each fit. The values in the remaining rows correspond to the parameters - for the supplied fit function. + for the supplied fit function, then the errors for each of the parameters. Failed fits have an rsquared of ``numpy.nan`` and parameters are set to the provided or constructed guess or ``numpy.nan`` - if no guess was provided or constructed. + if no guess was provided or constructed; errors are set to ``numpy.nan``. Raises ------ @@ -816,33 +853,44 @@ def image_fit(images, grid=None, function=gaussian2d, guess=None, plot=False): # Number of fit parameters the function accepts (minus 1 for xy). param_count = function.__code__.co_argcount - 1 - # Number of parameters to return (plus 1 for rsquared). - result_count = param_count + 1 + # Number of parameters to return: fitted parameters, errors, and plus 1 for rsquared. + result_count = 2 * param_count + 1 result = np.full((image_count, result_count), np.nan) # Construct guesses. - if guess is None: + if guess is None or guess is True: if function is gaussian2d: images_normalized = image_normalize(images, remove_field=True) centers = image_positions(images_normalized, grid=grid, normalize=False) variances = image_variances(images_normalized, centers=centers, grid=grid, normalize=False) - print(centers) - print(variances) + maxs = np.amax(images, axis=(1, 2)) mins = np.amin(images, axis=(1, 2)) guess = np.vstack(( centers, maxs - mins, mins, - np.sqrt(variances[0:2, :]), + np.sqrt(variances[:2, :]), variances[2, :] - )).transpose() + )).T else: - raise NotImplementedError("Default guess for function not implemented.") + message = f"Default guess for function {str(function)} not implemented." + if guess is True: + raise NotImplementedError(message) + else: + warnings.warn(message) # Fit and plot each image. for img_idx in range(image_count): img = images[img_idx, :, :].ravel() + grid_ravel_ = grid_ravel + + # Deal with nans. + undefined = np.isnan(img) + if np.any(undefined): + defined = np.logical_not(undefined) + img = img[defined] + grid_ravel_ = (grid_ravel[0][defined], grid_ravel[1][defined]) # Get guess. p0 = None if guess is None else guess[img_idx] @@ -850,9 +898,11 @@ def image_fit(images, grid=None, function=gaussian2d, guess=None, plot=False): # Attempt fit. fit_succeeded = True popt = None + perr = None try: - popt, _ = curve_fit(function, grid_ravel, img, ftol=1e-5, p0=p0,) + popt, pcov = curve_fit(function, grid_ravel_, img, ftol=1e-5, p0=p0,) + perr = np.sqrt(np.diag(pcov)) except RuntimeError: # The fit failed if scipy says so. fit_succeeded = False else: # The fit failed if any of the parameters aren't finite. @@ -860,15 +910,18 @@ def image_fit(images, grid=None, function=gaussian2d, guess=None, plot=False): fit_succeeded = False if fit_succeeded: # Calculate r2. - ss_res = np.sum(np.square(img - function(grid_ravel, *popt))) + ss_res = np.sum(np.square(img - function(grid_ravel_, *popt))) ss_tot = np.sum(np.square(img - np.mean(img))) r2 = 1 - (ss_res / ss_tot) else: # r2 is nan and the fit parameters are the guess or nan. popt = p0 if p0 is not None else np.full(param_count, np.nan) r2 = np.nan + perr = np.nan + # Populate results. result[img_idx, 0] = r2 - result[img_idx, 1:] = popt + result[img_idx, 1:(param_count+1)] = popt + result[img_idx, (param_count+1):] = perr # Plot. if plot: diff --git a/slmsuite/holography/analysis/fitfunctions.py b/slmsuite/holography/analysis/fitfunctions.py index 4323f8d..bea3ca0 100644 --- a/slmsuite/holography/analysis/fitfunctions.py +++ b/slmsuite/holography/analysis/fitfunctions.py @@ -4,6 +4,9 @@ import numpy as np + +# 1D + def linear(x, m, b): r""" For fitting a line. @@ -27,6 +30,31 @@ def linear(x, m, b): return m * x + b +def parabola(x, a, x0, y0): + r""" + For fitting a parabola. + + .. math:: y(x) = a(x - x_0)^2 + y_0 + + Parameters + ---------- + x : numpy.ndarray + Some independent variable. + a : float + Strength of the parabola. + x0 : float + :math:`x` position of the vertex. + y0 : float + :math:`y` position of the vertex. + + Returns + ------- + y : numpy.ndarray + Line evaluated at all ``x``. + """ + return a * np.square(x - x0) + y0 + + def hyperbola(z, w0, z0, zr): r""" For fitting a hyperbola. @@ -83,7 +111,7 @@ def lorentzian(x, x0, a, c, Q): r""" For fitting a resonance. There are many ways to describe a Lorentzian. Commonly, a full-width-half-maximum - definition is used. Here, with roots in photonic crystal cavities, we include a + definition is used. Here, with roots in photonic crystal cavities, we use a form defined using the quality factor :math:`Q` of the resonance. .. math:: y(x) = c + \frac{a}{1 + \left[\frac{x - x_0}{x_0/2Q}\right]^2}. @@ -181,6 +209,8 @@ def gaussian(x, x0, a, c, w): return c + a * np.exp(-.5 * np.square((x - x0) * (1/w))) +# 2D + def gaussian2d(xy, x0, y0, a, c, wx, wy, wxy=0): r""" For fitting a 2D Gaussian. diff --git a/slmsuite/holography/toolbox/__init__.py b/slmsuite/holography/toolbox/__init__.py index f62385e..d09402f 100644 --- a/slmsuite/holography/toolbox/__init__.py +++ b/slmsuite/holography/toolbox/__init__.py @@ -3,7 +3,7 @@ """ import numpy as np -from scipy.spatial.distance import chebyshev, euclidean +from scipy.spatial import distance from scipy.spatial import Voronoi, voronoi_plot_2d import cv2 import matplotlib.pyplot as plt @@ -28,16 +28,16 @@ CAMERA_UNITS = ["ij"] BLAZE_LABELS = { - "rad": (r"$\theta_x$ [rad]", r"$\theta_y$ [rad]"), + "rad": (r"$\theta_x$ [rad]", r"$\theta_y$ [rad]"), "mrad": (r"$\theta_x$ [mrad]", r"$\theta_y$ [mrad]"), - "deg": (r"$\theta_x$ [$^\circ$]", r"$\theta_y$ [$^\circ$]"), + "deg": (r"$\theta_x$ [$^\circ$]", r"$\theta_y$ [$^\circ$]"), "norm": (r"$k_x/k$", r"$k_y/k$"), - "kxy": (r"$k_x/k$", r"$k_y/k$"), - "knm": (r"$n$ [pix]", r"$m$ [pix]"), + "kxy": (r"$k_x/k$", r"$k_y/k$"), + "knm": (r"$k_n$ [pix]", r"$k_m$ [pix]"), "freq": (r"$f_x$ [1/pix]", r"$f_y$ [1/pix]"), "lpmm": (r"$k_x/2\pi$ [1/mm]", r"$k_y/2\pi$ [1/mm]"), - "ij": (r"Camera $x$ [pix]", r"Camera $y$ [pix]"), - "zernike": (r"$Z_1^1$ [Zernike rad]", r"$Z_1^{-1}$ [Zernike rad]"), + "ij": (r"Camera $x$ [pix]", r"Camera $y$ [pix]"), + "zernike": (r"$Z_1^1$ [Zernike rad]", r"$Z_1^{-1}$ [Zernike rad]"), } for k in LENGTH_FACTORS.keys(): u = LENGTH_LABELS[k] @@ -229,7 +229,7 @@ def convert_vector(vector, from_units="norm", to_units="norm", slm=None, shape=N cameraslm = None if from_units in CAMERA_UNITS or to_units in CAMERA_UNITS: - if cameraslm is None or cameraslm.fourier_calibration is None: + if cameraslm is None or cameraslm.calibrations["fourier"] is None: warnings.warn( f"CameraSLM must be passed to slm for conversion '{from_units}' to '{to_units}'" ) @@ -425,7 +425,7 @@ def convert_radius(radius, from_units="norm", to_units="norm", slm=None, shape=N def window_slice(window, shape=None, centered=False, circular=False): """ - Get the slices that describe the window's view into the larger array. + Parses the slices that describe the window's view into the larger array. Parameters ---------- @@ -476,14 +476,17 @@ def window_slice(window, shape=None, centered=False, circular=False): xc = xi + int((window[1] - 1) / 2) yc = yi + int((window[3] - 1) / 2) - rr_grid = (window[3] ** 2) * np.square(x_grid.astype(float) - xc) + ( - window[1] ** 2 - ) * np.square(y_grid.astype(float) - yc) + rr_grid = ( + (window[3] ** 2) * np.square(x_grid.astype(float) - xc) + + (window[1] ** 2) * np.square(y_grid.astype(float) - yc) + ) mask_grid = rr_grid <= (window[1] ** 2) * (window[3] ** 2) / 4.0 + # Pass things back through window_slice to crop the circle, should the user + # have given values that are out of bounds. return window_slice((y_grid[mask_grid], x_grid[mask_grid]), shape=shape) - else: # Otherwise, return square slices + else: # Otherwise, return square slices in the python style. slice_ = (slice(yi, yf), slice(xi, xf)) # (y_ind, x_ind) format elif len(window) == 2: @@ -714,6 +717,7 @@ def imprint( matrix : numpy.ndarray The data to imprint a ``function`` onto. window + Passed to :meth:`~slmsuite.holography.toolbox.window_slice()`. See :meth:`~slmsuite.holography.toolbox.window_slice()` for various options. function : function OR float A function in the style of :mod:`~slmsuite.holography.toolbox` helper functions, @@ -726,7 +730,7 @@ def imprint( Note ~~~~ - Also accepts floating point values, in which case this value is simply added. + Also accepts floating point values, in which case the value is simply added. grid : (array_like, array_like) OR :class:`~slmsuite.hardware.slms.slm.SLM` OR None Meshgrids of normalized :math:`\frac{x}{\lambda}` coordinates corresponding to SLM pixels, in ``(x_grid, y_grid)`` form. @@ -1080,15 +1084,15 @@ def fit_3pt(y0, y1, y2, N=None, x0=(0, 0), x1=(1, 0), x2=(0, 1), orientation_che return np.matmul(M, indices) + b -def smallest_distance(vectors, metric=chebyshev): +def smallest_distance(vectors, metric="chebyshev"): """ Returns the smallest distance between pairs of points under a given ``metric``. - Note - ~~~~ - An :math:`\mathcal{O}(N^2)` brute force approach is currently implemented. - Future work will involve an :math:`\mathcal{O}(N\log(N))` - divide and conquer algorithm. + Tip + ~~~ + This function supports a :math:`\mathcal{O}(N\log(N))` divide and conquer algorithm + and can handle large pointsets. + An :math:`\mathcal{O}(N^2)` brute force approach is implemented as a backup. Caution ~~~~~~~ @@ -1100,10 +1104,13 @@ def smallest_distance(vectors, metric=chebyshev): vectors : array_like Points to compare. Cleaned with :meth:`~slmsuite.holography.toolbox.format_2vectors()`. - metric : function + metric : str OR function Function to use to compare. - Defaults to :meth:`scipy.spatial.distance.chebyshev()`. - :meth:`scipy.spatial.distance.euclidean()` is also common. + Defaults to ``"chebyshev"`` which corresponds to + :meth:`scipy.spatial.distance.chebyshev()`. + The :math:`\mathcal{O}(N\log(N))` divide and conquer algorithm is only + compatible with string inputs. allowed by :meth:`scipy.spatial.distance.pdist`. + Function arguments will fallback to the brute force approach. Returns ------- @@ -1111,21 +1118,73 @@ def smallest_distance(vectors, metric=chebyshev): Minimum distance between any pair of points under the given metric. If fewer than two points are given, then ``np.inf`` is returned. """ + + def _divide_and_conquer_recursive(v, metric, axis=0, min_div=200): + # Expects sorted v. + N = v.shape[0] + + if N > min_div: + M = int(N/2) + + # Divide the problem recursively. + d1 = _divide_and_conquer_recursive(v[:M, :], metric, axis) + d2 = _divide_and_conquer_recursive(v[M:, :], metric, axis) + + # Conquer. + d = min(d1, d2) + + # Leave if we don't need to merge. + if (v[M, axis] - v[M+1, axis]) > d: + return d + + # Merge around average x0 between two sections. + x0 = (v[M, axis] + v[M+1, axis]) / 2 + mask = np.abs(v[:, axis] - x0) < d + subset = v[mask, :] + + return min(d, distance.pdist(subset, metric=metric).min()) + else: + # Use pdist as a fast low-level distance calculator. + return distance.pdist(v, metric=metric).min() + vectors = format_2vectors(vectors) N = vectors.shape[1] - minimum = np.inf - if N <= 1: - return minimum + return np.inf + + if isinstance(metric, str): # Divide and conquer. + if not metric in distance._METRIC_ALIAS: + raise RuntimeError("Distance metric '{metric}' not recognized by scipy.") + + axis = 0 + min_div = 200 + + # pdist needs transpose. + vectors = vectors.T - for x in range(N - 1): - for y in range(x + 1, N): - distance = metric(vectors[:, x], vectors[:, y]) - if distance < minimum: - minimum = distance + if N < 2*min_div: + return distance.pdist(vectors, metric=metric).min() + else: + centroid = np.max(vectors, axis=axis, keepdims=True) + + # Slightly inefficient use of cdist. + xorder = distance.cdist(vectors[:,[axis]], centroid[:,[axis]], metric=metric) + + I = np.argsort(np.squeeze(xorder)) + vsort = vectors[I, :] - return minimum + return _divide_and_conquer_recursive(vsort, metric, axis=axis, min_div=min_div) + else: # Fallback to brute force. + minimum = np.inf + + for x in range(N - 1): + for y in range(x + 1, N): + distance = metric(vectors[:, x], vectors[:, y]) + if distance < minimum: + minimum = distance + + return minimum def lloyds_algorithm(grid, vectors, iterations=10, plot=False): @@ -1355,7 +1414,7 @@ def transform_grid(grid, transform=None, shift=None, direction="fwd"): if not np.isscalar(transform): transform = np.squeeze(transform) if transform.shape != (2,2): - raise ValueError("Expected transform to be None, scalar or a 2x2 matrix.") + raise ValueError("Expected transform to be None, scalar, or a 2x2 matrix.") # Parse shift. if shift is None: @@ -1426,9 +1485,8 @@ def pad(matrix, shape): (shape[1] - matrix.shape[1]) / 2.0, ) - assert ( - deltashape[0] >= 0 and deltashape[1] >= 0 - ), "Shape {} is too large to pad to shape {}".format(tuple(matrix.shape), shape) + if not (deltashape[0] >= 0 and deltashape[1] >= 0): + raise ValueError(f"Shape {tuple(matrix.shape)} is too large to pad to shape {shape}") pad_b = int(np.floor(deltashape[0])) pad_t = int(np.ceil(deltashape[0])) @@ -1477,9 +1535,8 @@ def unpad(matrix, shape): deltashape = ((shape[0] - mshape[0]) / 2.0, (shape[1] - mshape[1]) / 2.0) - assert ( - deltashape[0] <= 0 and deltashape[1] <= 0 - ), "Shape {} is too small to unpad to shape {}".format(tuple(mshape), shape) + if not (deltashape[0] <= 0 and deltashape[1] <= 0): + raise ValueError(f"Shape {tuple(mshape)} is too small to unpad to shape {shape}") pad_b = int(np.floor(-deltashape[0])) pad_t = int(mshape[0] - np.ceil(-deltashape[0])) diff --git a/slmsuite/holography/toolbox/cuda.cu b/slmsuite/holography/toolbox/cuda.cu index dac1470..6eb69e9 100644 --- a/slmsuite/holography/toolbox/cuda.cu +++ b/slmsuite/holography/toolbox/cuda.cu @@ -433,22 +433,14 @@ extern "C" __global__ void update_weights_generic( int i = blockDim.x * blockIdx.x + threadIdx.x; if (i < N) { - float feedback = feedback_amp[i]; - float target = 0; + float feedback = feedback_amp[i] / feedback_norm; + float target = target_amp[i] - if (feedback_norm != 0) { - feedback /= feedback_norm; - target = target_amp[i] - } else { - target = 1 / sqrt(N); - } - - if (method != 4 && method != 5) { - // Handle feedback with a non-uniform target + if (method != 4 && method != 5) { // Multiplicative. if (target != 0) { feedback /= target; } else { - feedback = 1; // Do nothing. + feedback = 1; } if (method == 1 || method == 2) { // Leonardo, Kim @@ -456,18 +448,18 @@ extern "C" __global__ void update_weights_generic( } else if (method == 3) { // Nogrette feedback = 1 / (1 - feedback_factor * (1 - feedback)) } - } else { + } else { // Additive. if (method == 4) { // Wu feedback = exp(feedback_exponent * (target - feedback)) } else if (method == 5) { // tanh - feedback = feedback_factor * tanh(feedback_exponent * (target - feedback)) + feedback = 1 + feedback_factor * tanh(feedback_exponent * (target - feedback)) } } // Check nan, inf if (isinf(feedback) || isnan(feedback)) { feedback = 1 } - // Export the result to global memory. + // Export the result to global memory if changed. if (feedback != 1) { weight_amp[i] *= feedback; } diff --git a/slmsuite/holography/toolbox/phase.py b/slmsuite/holography/toolbox/phase.py index d14e322..09f2e51 100644 --- a/slmsuite/holography/toolbox/phase.py +++ b/slmsuite/holography/toolbox/phase.py @@ -22,12 +22,12 @@ # Basic gratings. -def blaze(grid, vector=(0, 0), offset=0): +def blaze(grid, vector=(0, 0)): r""" Returns a simple `blazed grating `_, a linear phase ramp, toward a given vector in :math:`k`-space. - .. math:: \phi(\vec{x}) = 2\pi \cdot \vec{k} \cdot \vec{x} + o + .. math:: \phi(\vec{x}) = 2\pi \cdot \vec{k} \cdot \vec{x} Parameters ---------- @@ -39,8 +39,6 @@ def blaze(grid, vector=(0, 0), offset=0): vector : (float, float) :math:`\vec{k}`. Blaze vector in normalized :math:`\frac{k_x}{k}` units. See :meth:`~slmsuite.holography.toolbox.convert_vector()` - offset : float - Phase offset for this blaze. Returns ------- @@ -59,25 +57,21 @@ def blaze(grid, vector=(0, 0), offset=0): else: result = (2 * np.pi * vector[0]) * x_grid + (2 * np.pi * vector[1]) * y_grid - # Add offset if provided. - if offset != 0: - result += offset - return result -def sinusoid(grid, vector=(0, 0), offset=0, amplitude=np.pi, outer_offset=np.pi): +def sinusoid(grid, vector=(0, 0), shift=0, a=np.pi, b=0): r""" Returns a simple `holographic grating `_, a sinusoidal grating, toward a given vector in :math:`k`-space. - .. math:: \phi(\vec{x}) = a \sin(\vec{k} \cdot \vec{x} + o) + b + .. math:: \phi(\vec{x}) = \frac{a-b}{2} [1 + \cos(2\pi \cdot \vec{k} \cdot \vec{x} + s)] + b Important --------- - Unlike a blazed grating :meth:`.blaze()`, half the power will be deflected toward - the mirror -1st order at :math:`-\vec{k}`, at best. + Unlike a blazed grating :meth:`.blaze()`, power will efficiently be deflected toward + the mirror -1st order at :math:`-\vec{k}`, by symmetry. Parameters ---------- @@ -89,13 +83,15 @@ def sinusoid(grid, vector=(0, 0), offset=0, amplitude=np.pi, outer_offset=np.pi) vector : (float, float) :math:`\vec{k}`. Blaze vector in normalized :math:`\frac{k_x}{k}` units. See :meth:`~slmsuite.holography.toolbox.convert_vector()` - offset : float - Grating phase offset. - amplitude : float - Amplitude of the sinusoid. - The 0th order will be minimized when this is equal to :math:`\pi`. - outer_offset : float - Phase offset for this grating. + shift : float + Radians to laterally shift the period of the grating by. + a : float + Value at one extreme of the sinusoid. + Ignoring crosstalk, + the 0th order will be minimized when ``|a-b|`` is equal to :math:`\pi`. + b : float + Value at the other extreme of the sinusoid. + Defaults to zero, in which case ``a`` is the amplitude. Returns ------- @@ -104,26 +100,26 @@ def sinusoid(grid, vector=(0, 0), offset=0, amplitude=np.pi, outer_offset=np.pi) """ if vector[0] == 0 and vector[1] == 0: (x_grid, _) = _process_grid(grid) - result = np.full_like(x_grid, amplitude * np.sin(offset)) + result = np.full_like(x_grid, (a-b)/2 * (1 + np.sin(shift))) else: - result = amplitude * np.sin(blaze(grid, vector, offset)) + result = (a-b)/2 * (1 + np.sin(blaze(grid, vector) + shift)) # Add offset if provided. - if outer_offset != 0: - result += outer_offset + if b != 0: + result += b return result -def binary(grid, vector=(0, 0), offset=0, amplitude=np.pi, outer_offset=0, duty_cycle=.5): +def binary(grid, vector=(0, 0), shift=0, a=np.pi, b=0, duty_cycle=.5): r""" Returns a simple binary grating toward a given vector in :math:`k`-space. .. math:: \phi(\vec{x}) = \left\{ \begin{array}{ll} - a+b, & ( - [2\pi \cdot \vec{k} \cdot \vec{x} + o] \,\,\,\,\text{mod}\,\,\,\, 2\pi + a, & ( + [2\pi \cdot \vec{k} \cdot \vec{x} + s] \,\,\,\,\text{mod}\,\,\,\, 2\pi ) < 2\pi*d \\ b, & \text{ otherwise}. \end{array} @@ -146,13 +142,13 @@ def binary(grid, vector=(0, 0), offset=0, amplitude=np.pi, outer_offset=0, duty_ vector : (float, float) :math:`\vec{k}`. Blaze vector in normalized :math:`\frac{k_x}{k}` units. See :meth:`~slmsuite.holography.toolbox.convert_vector()` - offset : float - Grating phase offset. - amplitude : float - Amplitude of the sinusoid. - The 0th order will be minimized when this is equal to :math:`\pi`. - outer_offset : float - Phase offset for this grating. + shift : float + Radians to laterally shift the period of the grating by. + a : float + Value at one extreme of the binary grating. + b : float + Value at the other extreme of the binary grating. + Defaults to zero, in which case ``a`` is the amplitude. duty_cycle : float Ratio of the period which is 'on'. @@ -165,10 +161,10 @@ def binary(grid, vector=(0, 0), offset=0, amplitude=np.pi, outer_offset=0, duty_ if vector[0] == 0 and vector[1] == 0: (x_grid, _) = _process_grid(grid) - phase = outer_offset - if offset != 0: - if np.mod(offset, 2*np.pi) < (2 * np.pi * duty_cycle): - phase += amplitude + phase = b + if shift != 0: + if np.mod(shift, 2*np.pi) < (2 * np.pi * duty_cycle): + phase = a result = np.full(x_grid.shape, phase) elif vector[0] != 0 and vector[1] != 0: pass # xor the next case. @@ -185,9 +181,9 @@ def binary(grid, vector=(0, 0), offset=0, amplitude=np.pi, outer_offset=0, duty_ # If we have not set result, then we have to use the slow np.mod option. if result is None: result = np.where( - np.mod(blaze(grid, vector, offset), 2*np.pi) < (2 * np.pi * duty_cycle), - offset, - offset + amplitude, + np.mod(blaze(grid, vector, shift), 2*np.pi) < (2 * np.pi * duty_cycle), + b, + a, ) return result @@ -758,7 +754,7 @@ def _plot_zernike_pyramid(grid, order, scale=1, **kwargs): **kwargs Passed to :meth:`.zernike()`. """ - order += 1 + order = int(order + 1) indices_ansi = np.arange((order * (order + 1)) // 2) indices_radial = zernike_convert_index(indices_ansi, from_index="ansi", to_index="radial") @@ -1216,21 +1212,26 @@ def _determine_source_radius(grid, w=None): The radius of the phase pattern in normalized :math:`\frac{x}{\lambda}` units. To produce perfect structured beams, this radius is equal to the radius of the gaussian profile of the source (ideally not clipped by the SLM). + If an SLM was passed as grid, retrieves the data from + :attr:`slmsuite.hardware.slms.slm.SLM.source` and + :meth:`slmsuite.hardware.slms.slm.SLM.fit_source_amplitude()`. If ``w`` is left as ``None``, ``w`` is set to a quarter of the smallest normalized screen dimension. Returns ------- - float + w : float Determined radius. In normalized units. """ - (x_grid, y_grid) = _process_grid(grid) + if w is not None: + return w - # TODO + if hasattr(grid, "slm"): + grid = grid.slm + if hasattr(grid, "get_source_radius"): + return grid.get_source_radius() - if w is None: - return np.min([np.amax(x_grid), np.amax(y_grid)]) / 4 - else: - return w + (x_grid, y_grid) = _process_grid(grid) + return np.min([np.amax(x_grid), np.amax(y_grid)]) / 4 def laguerre_gaussian(grid, l, p, w=None):