[feat] added peak simulation

Author: yxkdejong

src/odemis/driver/simcam.py

@@ -34,9 +34,11 @@
 from odemis import dataio, model, util
 from odemis.model import oneway
 from odemis.util.synthetic import ParabolicMirrorRayTracer
+from odemis.util.synthetic import simulate_peak
 
 ERROR_STATE_FILE = "simcam-hw.error"
-
+GOFFSET_TO_PIXEL = 0.25  # Conversion factor for grating offset to image pixels.
+PEAK_WIDTH = 2.5  # Width of the simulated spectrograph peak in pixels (before binning).
 
 class Camera(model.DigitalCamera):
     '''
@@ -179,6 +181,22 @@ def clip_max_res(img_res):
             self._img_simulator = None
             self._mirror = None
 
+        try:
+            spectrograph = dependencies["spectrograph"]
+            if not (
+                    isinstance(spectrograph, model.ComponentBase)
+                    and hasattr(spectrograph, "axes")
+                    and isinstance(spectrograph.axes, dict)
+                    and "goffset" in spectrograph.axes):
+                raise ValueError("spectrograph %s must have a 'goffset' attribute" % spectrograph)
+
+            self._spectrograph = spectrograph
+            logging.debug("Will simulate spectral peaks using spectrograph %s", spectrograph.name)
+
+        except (TypeError, KeyError):
+            logging.info("Will not simulate spectrograph peaks")
+            self._spectrograph = None
+
         # Simple implementation of the flow: we keep generating images and if
         # there are subscribers, they'll receive it.
         self.data = SimpleDataFlow(self)
@@ -374,6 +392,11 @@ def _simulate(self):
         center = self._img_res[0] / 2, self._img_res[1] / 2
         pixel_size = self._metadata.get(model.MD_PIXEL_SIZE, self.pixelSize.value)
 
+        # Extra translation to simulate stage movement
+        pos = self._metadata.get(model.MD_POS, (0, 0))
+        pxs = [p / b for p, b in zip(pixel_size, self.binning.value)]
+        stage_shift = pos[0] / pxs[0], -pos[1] / pxs[1]  # Y goes opposite
+
         if self._mirror:
             eff_pixel_size = [p * b for p, b in zip(pixel_size, binning)]
             self._img_simulator.resolution.value = res
@@ -382,11 +405,6 @@ def _simulate(self):
             self._img = model.DataArray(sim_img, self._img.metadata)
             sim_img = self._img.copy()
         else:
-            # Extra translation to simulate stage movement
-            pos = self._metadata.get(model.MD_POS, (0, 0))
-            pxs = [p / b for p, b in zip(pixel_size, self.binning.value)]
-            stage_shift = pos[0] / pxs[0], -pos[1] / pxs[1]  # Y goes opposite
-
             # First and last index (eg, 0 -> 255)
             ltrb = [center[0] + trans[0] + stage_shift[0] - (res[0] / 2) * binning[0],
                     center[1] + trans[1] + stage_shift[1] - (res[1] / 2) * binning[1],
@@ -413,6 +431,31 @@ def _simulate(self):
                      [int(round(ltrb[1] + i * binning[1])) for i in range(res[1])])
             sim_img = self._img[numpy.ix_(coord[1], coord[0])]  # copy
 
+        # spectrograph peak simulation
+        if self._spectrograph:
+            current_offset = self._spectrograph.position.value["goffset"]
+
+            ccd_center_x = self._img_res[0]/2.0  # find the x-coordinate of the center of the ccd
+            x0_px = ccd_center_x + current_offset * GOFFSET_TO_PIXEL
+            roi_left = center[0] + trans[0] + stage_shift[0] - (res[0] / 2) * binning[0]
+
+            bin_x = binning[0]  # binning factor along x-axis
+            peak_center_binned = (x0_px - roi_left) / bin_x  # express the peak position in the ROI's coordinate system
+
+            logging.info("DEBUG: x0_px=%s, ltrb0=%s, result=%s",
+                         x0_px, roi_left, peak_center_binned
+                         )
+
+            width_binned = PEAK_WIDTH / bin_x
+
+            peak = simulate_peak(amplitude=20000, x0=peak_center_binned, width=width_binned,
+                                shape=sim_img.shape, dtype=sim_img.dtype)
+
+            # set all values in sim_img to the minimal sim_img value
+            min_val = sim_img.min()
+            sim_img[...] = min_val
+            sim_img += peak
+
         # Add some noise
         mx = self._img.max()
         sim_img += numpy.random.randint(0, max(mx // 100, 10), sim_img.shape, dtype=sim_img.dtype)

src/odemis/driver/test/simcam_test.py

@@ -457,5 +457,79 @@ def test_acquire_ar_pol(self):
         testing.assert_array_not_equal(im_horizontal, im_vertical)
 
 
+class TestSimCamSpectrograph(unittest.TestCase):
+
+    @classmethod
+    def setUpClass(cls):
+        class MockSpectrograph(model.Component):
+            def __init__(self, name="mock_spectrograph"):
+                super().__init__(name)
+                self.axes = {"goffset": None}
+                self.position = type('obj', (object,), {'value': {"goffset": 0.0}})()
+
+        cls.focus = simulated.Stage(**KWARGS_FOCUS)
+        cls.mock_spectrograph = MockSpectrograph()
+
+        cls.camera = CLASS(dependencies={"focus": cls.focus, "spectrograph": cls.mock_spectrograph},
+                               **KWARGS_MOVE)
+
+    @classmethod
+    def tearDownClass(cls):
+        cls.camera.terminate()
+
+    def setUp(self):
+        self.camera.binning.value = (1, 1)
+        self.camera.translation.value = (0, 0)
+        self.camera.resolution.value = self.camera.resolution.range[1]
+        self.camera.exposureTime.value = 0.05
+        self.mock_spectrograph.position.value["goffset"] = 0.0
+        time.sleep(0.1)
+
+    # helper function to find the peak
+    def _find_peak(self, image):
+        data = image.astype(float)
+        profile = data.mean(axis=0)
+
+        p_min, p_max = profile.min(), profile.max()
+        threshold = p_min + (p_max - p_min) * 0.8
+
+        mask = profile > threshold
+        indices = numpy.arange(len(profile))
+
+        if numpy.sum(mask) == 0:
+            return float(profile.argmax(()), p_max, profile.mean())
+
+        # find the peak
+        weighted_average = numpy.sum(indices[mask] * profile[mask]) / numpy.sum(profile[mask])
+
+        return weighted_average, p_max, profile.mean()
+
+    def test_spectrograph_with_goffset(self):
+        move = 300
+        expected_ratio = 0.25
+
+        self.mock_spectrograph.position.value["goffset"] = 0.0
+        image_0 = self.camera.data.get()
+        date_0 = image_0.metadata[model.MD_ACQ_DATE]
+
+        # change offset
+        self.mock_spectrograph.position.value["goffset"] = move
+
+        image_new = None
+        for _ in range(10):
+            image_new = self.camera.data.get()
+            if image_new.metadata[model.MD_ACQ_DATE] > date_0:
+                break
+            time.sleep(0.05)
+
+        x0, _, _ = self._find_peak(image_0)
+        x1, _, _ = self._find_peak(image_new)
+
+        shift = x1 - x0
+        print(f"DEBUG: Moved {move}, Peak went from {x0} to {x1}. Shift: {shift}")
+
+        expected_shift = move * expected_ratio
+        self.assertAlmostEqual(shift, expected_shift, delta=0.5)
+
 if __name__ == '__main__':
     unittest.main()

src/odemis/util/synthetic.py

@@ -128,6 +128,37 @@ def psf_gaussian(
     numpy.rint(UINT16_MAX * out, out=out)
     return out.astype(numpy.uint16)
 
+def simulate_peak(amplitude, x0, width, shape, dtype=numpy.uint16) -> numpy.ndarray:
+
+    """
+    Simulate a 1D Gaussian peak across the x-dimension of an image, with optional 2D image input for direct use when y > 1.
+    :param amplitude: the maximum intensity of the peak
+    :param x0: the center position of the peak along the x-axis
+    :param width: the standard deviation (width) of the Gaussian peak
+    :param shape: the shape of the output image (can be an int for 1D or a tuple for 2D)
+    :param dtype: the data type of the output array (default is numpy.uint16)
+    :return: a numpy array containing the simulated peak, either as a 1D Gaussian or directly from the provided image
+    """
+
+    # Standardize shape to a tuple, so we can index it
+    if isinstance(shape, (int, numpy.integer)):
+        shape = (shape,)
+
+    # Generate the 1D Gaussian peak across the x-dimension
+    x = numpy.arange(shape[-1])
+    intensity = amplitude*(numpy.exp(-0.5 * ((x - x0) / width)**2))
+
+    # Clip based on the actual dtype provided
+    info = numpy.iinfo(dtype)
+    intensity_clipped = numpy.clip(intensity, info.min, info.max)
+    peak_1d = intensity_clipped.astype(dtype)
+
+    if len(shape) == 1:
+        return peak_1d
+    else:
+        # Duplicated along the y-dimension (shape[0])
+        return numpy.tile(peak_1d, (shape[0], 1))
+
 
 class ParabolicMirrorRayTracer:
     """

src/odemis/util/test/synthetic_test.py

@@ -23,7 +23,7 @@
 
 import numpy
 
-from odemis.util.synthetic import ParabolicMirrorRayTracer
+from odemis.util.synthetic import ParabolicMirrorRayTracer, simulate_peak
 
 
 class TestParabolicMirrorRayTracer(unittest.TestCase):
@@ -125,6 +125,84 @@ def test_constructor_raises_value_error_on_missing_keys(self):
         with self.assertRaises(ValueError):
             ParabolicMirrorRayTracer(good_pos=bad_pos_missing_all)
 
+class TestPeakSimulation(unittest.TestCase):
+
+    def test_simulate_peak_1d(self):
+        """
+        Verify simulate_peak generates a correct 1D Gaussian peak.
+        """
+
+        amplitude = 1000
+        x0 = 50
+        width = 5
+        shape = 100
+
+        peak = simulate_peak(amplitude, x0, width, shape)
+
+        self.assertEqual(peak.shape, (shape,))
+        self.assertEqual(peak.dtype, numpy.uint16)
+
+        # peak maximum should occur near x0
+        peak_idx = numpy.argmax(peak)
+        self.assertAlmostEqual(peak_idx, x0, delta=1)
+
+        # peak value should be close to amplitude (clipped if necessary)
+        self.assertAlmostEqual(peak[peak_idx], amplitude, delta=5)
+
+    def test_simulate_peak_2d(self):
+        """
+        Verify simulate_peak correctly expands the 1D peak to 2D.
+        """
+
+        amplitude = 500
+        x0 = 40
+        width = 4
+        shape = (20, 100)
+
+        peak = simulate_peak(amplitude, x0, width, shape)
+
+        self.assertEqual(peak.shape, shape)
+
+        # every row should be identical
+        for row in peak:
+            self.assertTrue(numpy.array_equal(row, peak[0]))
+
+        # verify peak location
+        peak_idx = numpy.argmax(peak[0])
+        self.assertAlmostEqual(peak_idx, x0, delta=1)
+
+    def test_simulate_peak_dtype_clipping(self):
+        """
+        Verify values are clipped to dtype limits.
+        """
+
+        amplitude = 100000  # larger than uint16 max
+        x0 = 30
+        width = 3
+        shape = 80
+
+        peak = simulate_peak(amplitude, x0, width, shape)
+
+        dtype_max = numpy.iinfo(numpy.uint16).max
+
+        self.assertEqual(peak.max(), dtype_max)
+
+    def test_simulate_peak_dtype(self):
+        """
+        Verify simulate_peak respects the dtype argument.
+        """
+
+        amplitude = 200
+        x0 = 25
+        width = 3
+        shape = 60
+
+        peak = simulate_peak(amplitude, x0, width, shape, dtype=numpy.uint8)
+
+        self.assertEqual(peak.dtype, numpy.uint8)
+
+        dtype_max = numpy.iinfo(numpy.uint8).max
+        self.assertLessEqual(peak.max(), dtype_max)
 
 if __name__ == "__main__":
     unittest.main()