Rewrite for the different backends happening here

photonics/experiments/backends.py

@@ -0,0 +1,85 @@
+# defines "photonic backends"
+# i.e. DM-WFS
+# the children yearn for the Julia type tree
+# actually why don't I just port this whole thing to Julia? DAO has an interface for it...
+# because it's too much work to rewrite everything
+
+from abc import ABC
+
+class PhotonicBackend(ABC):
+    """
+    - get intensities off an image given that we've got centroids and a radius
+    - measure darks and PL flats and save them appropriately
+    - take an interaction matrix and invert it
+    - make linearity curves
+    - do a "pseudo-closed-loop" and a real closed loop
+    """
+    def set_centroids(self, centroids, image, save=True):
+        """
+        Sets the centroids of each PL port based on 
+        - an initial guess of each position, in `centroids`
+        - a reference image, in `image`.
+        """
+        self.xc, self.yc = ports_in_radial_order(refine_centroids(centroids))
+        self.yi, self.xi = np.indices(image.shape)
+        self.masks = [(self.xi - xc) ** 2 + (self.yi - yc) ** 2 <= self.spot_radius_px ** 2 for (xc, yc) in zip(self.xc, self.yc)]
+        self.centroids_dt = datetime_now()
+        if save:
+            new_centroids = np.zeros_like(centroids)
+            new_centroids[:,0] = self.xc
+            new_centroids[:,1] = self.yc
+            with h5py.File(self.filepath(f"centroids_{self.centroids_dt}", ext="hdf5"), "w") as f:
+                c = f.create_dataset("centroids", data=new_centroids)
+                c.attrs["spot_radius_px"] = self.spot_radius_px
+
+    def get_intensities(self, img, exclude=[]):
+        intensities = np.zeros(self.Nports - len(exclude))
+        j = 0
+        for i in range(self.Nports - len(exclude)):
+            while j in exclude:
+                j += 1
+            intensities[i] = np.sum(img[self.masks[i]])
+            j += 1
+
+        return normalize(intensities)
+
+    def reconstruct_image(self, img, intensities):
+        recon_image = np.zeros_like(img)
+        for (i, intensity) in enumerate(intensities):
+            recon_image[self.masks[i]] = intensity
+
+        return recon_image
+
+    @property
+    def directory(self):
+        return path.join(DATA_PATH, self.subdir)
+
+    def filepath(self, fname, ext=None):
+        """
+        The path we want to save/load data from or to.
+        """
+        if ext is None:
+            ext = self.ext
+        return path.join(DATA_PATH, self.subdir, fname + "." + ext)
+
+    def measure_dark(self, direct=True):
+        if direct:
+            input("Taking a dark frame, remove the light source!")
+        darks = []
+        for _ in trange(self.nframes):
+            darks.append(self.im.get_data(check=True).astype(float))
+            sleep(self.exp_ms / 1000)
+
+        self.dark = np.mean(darks, axis=0)
+        self.save(f"dark_exptime_ms_{self.exp_ms}_gain_{self.gain}", self.dark)
+
+    def measure_pl_flat(self):
+        self.send_zeros(verbose=True)
+        self.pl_flat = self.get_intensities(self.get_image())
+        self.save(f"pl_flat_{datetime_now()}", self.pl_flat)
+
+class SimBackend(PhotonicBackend):
+    pass
+
+class muirSEALBackend(PhotonicBackend):
+    pass

photonics/experiments/lantern_reader.py

@@ -113,29 +113,6 @@ def plot_ports(self, save=False):
             plt.savefig(self.filepath(f"port_mask_{datetime_now()}", ext="png"))
         plt.show()
 
-    def ports_in_radial_order(self, points):
-        xnew, ynew = np.zeros_like(points[:,0]), np.zeros_like(points[:,1])
-        prev_idx = 0
-        radial_shell = np.zeros(len(xnew))
-        k = 0
-        while len(points) > 0: # should run three times for a 19-port lantern
-            if len(points) == 1:
-                v = np.array([0])
-            else:
-                hull = ConvexHull(points)
-                v = hull.vertices
-            nhull = len(v)
-            xtemp, ytemp = points[v][:,0], points[v][:,1]
-            sortperm = np.argsort(angles_relative_to_center(xtemp, ytemp))[::-1]
-            xnew[prev_idx:prev_idx+nhull] = xtemp[sortperm]
-            ynew[prev_idx:prev_idx+nhull] = ytemp[sortperm]
-            radial_shell[prev_idx:prev_idx+nhull] = k
-            k += 1
-            prev_idx += nhull
-            points = np.delete(points, v, axis=0)
-
-        return np.flip(xnew), np.flip(ynew), np.flip((k - 1 - radial_shell) / (k - 1))
-
     def bounding_box(self, pad=0):
         return (floor(min(self.xc)-pad), ceil(max(self.xc)+pad), floor(min(self.yc)-pad), ceil(max(self.yc)+pad))
 

photonics/experiments/seal.py

@@ -129,7 +129,7 @@ def monitor_stability(plwfs, cam, n=100):
 
     return measurements
 
-def (plwfs):
+def snr(plwfs):
     reader = plwfs.reader
     img = plwfs.cam.get()
     rp = 2 * (reader.fwhm + 1)

photonics/experiments/shane.py

@@ -8,7 +8,8 @@
 from tqdm import tqdm, trange
 import warnings
 import h5py
-from ..utils import date_now, datetime_now, datetime_ms_now, time_ms_now, rms, DATA_PATH, zernike_names
+
+from ..utils import date_now, datetime_now, datetime_ms_now, time_ms_now, rms, DATA_PATH, zernike_names, center_of_mass
 
 def normalize(x):
     return x / np.sum(x)
@@ -64,9 +65,18 @@ def setup_paramiko(self):
 
     def destroy_paramiko(self):
         self.client.close()
-
-    def set_centroids(self, centroids, save=True):
-        self.centroids = centroids
+
+    def refine_centroids(self, centroids, image, cutout_size=12):
+        new_centroids = np.zeros_like(centroids)
+        for i in range(self.Nports):
+            xl, xu = int(centroids[i,0]) - cutout_size, int(centroids[i,0]) + cutout_size + 1
+            yl, yu = int(centroids[i,1]) - cutout_size, int(centroids[i,1]) + cutout_size + 1
+            new_centroids[i] = center_of_mass(image[xl:xu, yl:yu], xg[xl:xu, yl:yu], yg[xl:xu, yl:yu])
+
+        return new_centroids
+
+    def set_centroids(self, centroids, image, save=True):
+        self.centroids = self.refine_centroids(centroids)
         self.xc, self.yc = centroids[:,0], centroids[:,1]
         self.masks = [(self.xi - xc) ** 2 + (self.yi - yc) ** 2 <= self.spot_radius_px ** 2 for (xc, yc) in self.centroids]
         self.centroids_dt = datetime_now()

photonics/utils.py

@@ -4,6 +4,7 @@
 from os import path
 from copy import copy
 from hcipy import Field, imshow_field
+from scipy.spatial import ConvexHull
 
 PROJECT_ROOT = path.dirname(path.dirname(path.abspath(__file__)))
 DATA_PATH = path.join(PROJECT_ROOT, "data")
@@ -71,3 +72,42 @@ def norm(a, b):
 
 def corr(a, b):
     return np.abs(norm(a, b)) / np.sqrt(norm(a, a) * norm(b, b))
+
+def center_of_mass(img, indices_x, indices_y):
+    # Returns the center of mass (intensity-weighted sum) in the x and y direction of the image.
+    xc = np.sum(img * indices_x) / np.sum(img)
+    yc = np.sum(img * indices_y) / np.sum(img)
+    return xc, yc
+
+
+def ports_in_radial_order(points):
+    xnew, ynew = np.zeros_like(points[:,0]), np.zeros_like(points[:,1])
+    prev_idx = 0
+    radial_shell = np.zeros(len(xnew))
+    k = 0
+    while len(points) > 0: # should run three times for a 19-port lantern
+        if len(points) == 1:
+            v = np.array([0])
+        else:
+            hull = ConvexHull(points)
+            v = hull.vertices
+        nhull = len(v)
+        xtemp, ytemp = points[v][:,0], points[v][:,1]
+        sortperm = np.argsort(angles_relative_to_center(xtemp, ytemp))[::-1]
+        xnew[prev_idx:prev_idx+nhull] = xtemp[sortperm]
+        ynew[prev_idx:prev_idx+nhull] = ytemp[sortperm]
+        radial_shell[prev_idx:prev_idx+nhull] = k
+        k += 1
+        prev_idx += nhull
+        points = np.delete(points, v, axis=0)
+
+    return np.flip(xnew), np.flip(ynew)
+
+def refine_centroids(centroids, image, cutout_size=12):
+    new_centroids = np.zeros_like(centroids)
+    for i in range(len(centroids)):
+        xl, xu = int(centroids[i,0]) - cutout_size, int(centroids[i,0]) + cutout_size + 1
+        yl, yu = int(centroids[i,1]) - cutout_size, int(centroids[i,1]) + cutout_size + 1
+        new_centroids[i] = center_of_mass(image[xl:xu, yl:yu], xg[xl:xu, yl:yu], yg[xl:xu, yl:yu])
+
+    return new_centroids

scripts_py/analysis/better_centroiding.py

@@ -0,0 +1,69 @@
+# %%
+import numpy as np
+from matplotlib import pyplot as plt
+# %%
+# it's hard to test if this is working as is...let's make up a test image
+
+def hexagon_pattern(nrings, core_offset):
+    nports = 1 + 3 * nrings * (nrings - 1)
+    port_positions = np.zeros((nports, 2))
+    nports_so_far = 0
+    for i in range(nrings):
+        nports_per_ring = max(1,6*i)
+        theta = 0
+        current_position = i * core_offset * np.array([np.cos(theta), np.sin(theta)])
+        next_position = i * core_offset * np.array([np.cos(theta + np.pi / 3), np.sin(theta + np.pi/3)])
+        for j in range(nports_per_ring):
+            if i > 0 and j % i == 0:
+                theta += np.pi / 3
+                current_position = next_position
+                next_position = i * core_offset * np.array([np.cos(theta + np.pi / 3), np.sin(theta + np.pi / 3)])
+            cvx_coeff = 0 if i == 0 else (j % i) / i
+            port_positions[nports_so_far,:] = (1 - cvx_coeff) * current_position + cvx_coeff * next_position
+            nports_so_far += 1
+    return port_positions
+# %%
+centroids = 300 + hexagon_pattern(3, 100)
+jitter = np.random.normal(0, 5, size=centroids.shape)
+centroids += jitter
+radii = np.random.uniform(4, 6, size=19)
+
+reference_image = np.zeros((520,656))
+xg, yg = np.indices((520,656))
+
+for ((xc, yc), r) in zip(centroids, radii):
+    reference_image += np.exp(-((xg - xc)**2 + (yg - yc)**2)/r**2)
+
+reference_image *= 4032 / np.max(reference_image)
+reference_image = np.int64(reference_image)
+reference_image += np.random.poisson(reference_image)
+
+plt.imshow(reference_image, cmap='gray')
+plt.show()
+# %%
+# now let's say I misidentified the centroids
+identified_centroids = np.random.normal(centroids, 3)
+plt.imshow(reference_image, cmap='gray')
+plt.scatter(identified_centroids[:,1], identified_centroids[:,0])
+plt.show()
+# %%
+def center_of_mass(img, indices_x, indices_y):
+    # Returns the center of mass (intensity-weighted sum) in the x and y direction of the image.
+    xc = np.sum(img * indices_x) / np.sum(img)
+    yc = np.sum(img * indices_y) / np.sum(img)
+    return xc, yc
+
+# %%
+cutout_size = 12
+new_centroids = np.zeros_like(centroids)
+for i in range(19):
+    xl, xu = int(identified_centroids[i,0]) - cutout_size, int(identified_centroids[i,0]) + cutout_size + 1
+    yl, yu = int(identified_centroids[i,1]) - cutout_size, int(identified_centroids[i,1]) + cutout_size + 1
+    new_centroids[i] = center_of_mass(reference_image[xl:xu, yl:yu], xg[xl:xu, yl:yu], yg[xl:xu, yl:yu])
+
+plt.imshow(reference_image, cmap='gray')
+plt.scatter(new_centroids[:,1], new_centroids[:,0], s=1)
+# %%
+print(float(np.std(identified_centroids - centroids)))
+print(float(np.std(new_centroids - centroids)))
+# %%

scripts_py/analysis/pl_241119_findspiral.py

@@ -0,0 +1,55 @@
+# %%
+import numpy as np
+from matplotlib import pyplot as plt
+import h5py
+from tqdm import trange
+from photonics.experiments.shane import ShaneLantern
+# %%
+with h5py.File("../../data/pl_241119/centroids_2024-11-19T18.19.22.185917.hdf5") as f:
+    centroids = np.array(f["centroids"])
+# %%
+lant = ShaneLantern()
+lant.set_centroids(centroids)
+
+# %%
+spiral_search_runs = [
+    # "stability_2024-11-19T18.35.31.751850.hdf5",
+    """    "stability_2024-11-19T18.49.20.552640.hdf5",
+    "stability_2024-11-19T18.52.19.198809.hdf5",
+    "stability_2024-11-19T18.52.28.898958.hdf5",
+    "stability_2024-11-19T18.53.36.906790.hdf5",
+    "stability_2024-11-19T18.53.46.656923.hdf5",
+    "stability_2024-11-19T18.54.44.361814.hdf5",
+    "stability_2024-11-19T18.57.36.228830.hdf5",
+    "stability_2024-11-19T19.02.49.030057.hdf5",
+    "stability_2024-11-19T19.06.03.842302.hdf5",
+    "stability_2024-11-19T19.08.08.360465.hdf5",
+    "stability_2024-11-19T19.13.02.699090.hdf5","""
+# %%
+spiral_search_runs = [
+    "stability_2024-11-19T21.08.09.370403.hdf5"
+]
+# %%
+with h5py.File(f"../../data/pl_241119/stability_2024-11-19T18.35.31.751850.hdf5") as f:
+    ref_image = np.array(f["pl_images"][0]) / 2
+# %%
+for x in spiral_search_runs:
+    with h5py.File(f"../../data/pl_241119/{x}") as f:
+        exp_ms = f["pl_intensities"].attrs["exp_ms"]
+        ref_image_conversion = exp_ms / 500
+        images = np.array(f["pl_images"])
+        import matplotlib.animation as animation
+
+        fig, ax = plt.subplots()
+        im = ax.imshow(ref_image * ref_image_conversion, cmap='viridis')
+
+        def update(frame):
+            im.set_array((images[frame, :, :] - ref_image * ref_image_conversion))
+            # (np.sum(lant.masks, axis=0)
+            return [im]
+
+        ani = animation.FuncAnimation(fig, update, frames=len(images), blit=True)
+        ani.save(f'animation_{x}.mp4', writer='ffmpeg')
+        plt.show()
+        plt.plot(np.array([lant.get_intensities(img - ref_image * ref_image_conversion) for img in images]))
+# %%