Fix sign and normalization in hybrid mode

deeptrack/optical/optics.py

@@ -890,7 +890,7 @@ def _pupil_numpy(
             * refractive_index_medium
             / wavelength
             * voxel_size[2]
-            * np.sqrt(1 - (NA / refractive_index_medium) ** 2 * RHO)
+            * (np.sqrt(1 - (NA / refractive_index_medium) ** 2 * RHO) - 1)
         )
 
         z_shift[z_shift.imag != 0] = 0
@@ -1088,7 +1088,7 @@ def _pupil_torch(
             inside,
             torch.zeros_like(inside),
         )
-        z_shift = (k0 * vz) * torch.sqrt(inside).to(complex_dtype)
+        z_shift = (k0 * vz) * (torch.sqrt(inside) - 1.0).to(complex_dtype)
 
         # defocus reshape (-1,1,1)
         defocus_t = torch.as_tensor(

deeptrack/optical/scatterers.py

@@ -2130,8 +2130,10 @@ def _solve_geometric(
             L, illumination_angle_field, coefficients
         )
 
+        # From Bohren & Huffman eq 4.75 E_s = E_0 * exp(i k R) / (-i k R) * S
+        # and 1 / (-i) = + i
         scattered_values = (
-            -1j
+            1j
             / (k * R3_field)
             * xp.exp(1j * k * R3_field)
             * (S2 * S2_coef + S1 * S1_coef)
@@ -2355,7 +2357,18 @@ def _solve_hybrid(
         )
         S1, S2 = self._mie_scattering(L, illum_valid, coefficients)
 
-        scattered_values = (S2 * S2_coef + S1 * S1_coef) / amp_factor
+        # Angular spectrum of the scattered field. A far field
+        # exp(i k R) / (-i k R) * S corresponds to plane-wave amplitudes
+        # A(kx, ky) = -S / (2 pi k^2 cos(theta)) (stationary phase / Weyl
+        # expansion). The DFT that follows needs A * (2 pi / dx)^2, giving
+        # -2 pi / ((k dx)^2 cos(theta)) * S relative to the unit illumination.
+        cos_valid = cos_theta_field[pupil_mask]
+        scattered_values = (
+            -2
+            * np.pi
+            / (k**2 * voxel_size[0] * voxel_size[1] * cos_valid)
+            * (S2 * S2_coef + S1 * S1_coef)
+        ) / amp_factor
 
         if TORCH_AVAILABLE and torch.is_tensor(arr):
             flat_values = torch.zeros(

tests/test_dlcc.py

@@ -72,9 +72,9 @@ def test_3_B(self):
 
         # First resolve
         expected_first = np.array(
-            [[[0.84823867], [0.84193078], [0.85179172]],
-            [[0.80131109], [0.79478238], [0.80499449]],
-            [[0.76094944], [0.75431347], [0.76469805]]]
+            [[[1.16725495], [1.17511633], [1.16286531]],
+            [[1.22730774], [1.23611139], [1.22238765]],
+            [[1.28311394], [1.29273146], [1.27773575]]]
         )
 
         with self._suppress_expected_optics_warnings():
@@ -88,9 +88,9 @@ def test_3_B(self):
 
         # Change after update
         expected_second = np.array(
-            [[[1.02220767], [0.98545219], [0.95156275]],
-            [[0.98707199], [0.94484236], [0.90718955]],
-            [[0.94221739], [0.89565235], [0.85521306]]]
+            [[[0.97935283], [1.01489583], [1.04966697]],
+            [[1.01348304], [1.05672786], [1.09800049]],
+            [[1.05950353], [1.11092199], [1.15915346]]]
         )
 
         with self._suppress_expected_optics_warnings():

tests/test_optics.py

@@ -112,6 +112,45 @@ def test_Brightfield(self):
         self.assertIsInstance(output_image, self.array_type)
         self.assertEqual(output_image.shape, (64, 64, 1))
 
+    def test_Brightfield_illumination_phase_independent_of_volume(self):
+        # The illumination must reach the camera as 1, wherever the sample
+        # volume sits in z. Otherwise its phase relative to fields added at
+        # the focal plane (e.g. MieSphere) depends on unrelated objects.
+        microscope = optics.Brightfield(
+            NA=1.0,
+            wavelength=633e-9,
+            resolution=100e-9,
+            magnification=1,
+            refractive_index_medium=1.33,
+            output_region=(0, 0, 64, 64),
+            return_field=True,
+        )
+        mie = MieSphere(
+            radius=0.3e-6,
+            refractive_index=1.58,
+            position=(32, 32),
+            position_unit="pixel",
+            z=0,
+            mode="hybrid",
+        )
+        far_sphere = Sphere(
+            radius=0.2e-6,
+            refractive_index=1.58,
+            position=(8, 8),
+            position_unit="pixel",
+            z=20,
+        )
+
+        with self._suppress_expected_optics_warnings():
+            sphere_only = np.asarray(microscope(far_sphere).resolve())[..., 0]
+            mie_only = np.asarray(microscope(mie).resolve())[..., 0]
+            both = np.asarray(microscope(mie >> far_sphere).resolve())[..., 0]
+
+        self.assertAlmostEqual(abs(sphere_only[56, 56] - 1), 0, delta=0.02)
+        window = (slice(22, 42), slice(22, 42))
+        change = np.abs(np.abs(both[window]) ** 2 - np.abs(mie_only[window]) ** 2)
+        self.assertLess(change.max(), 0.1)
+
     def test_Holography(self):
         microscope = optics.Holography(
             NA=0.7,

tests/test_scatterers.py

@@ -343,6 +343,58 @@ def test_MieSphere_modes(self):
         self.assertGreater(ratio, 1e-2)
         self.assertLess(ratio, 1e2)
 
+    def test_MieSphere_weak_sphere_forward_amplitude(self):
+        # A small, weakly scattering sphere is a thin phase object: summed
+        # over the image plane, its scattered field must equal
+        # +i k0 (n_p - n_m) V (Born limit, exp(-i w t)). The result must not
+        # depend on the simulation pixel size, and its sign means that a
+        # sphere denser than the medium delays the light.
+        wavelength, n_m, n_p, radius = 532e-9, 1.33, 1.40, 30e-9
+        expected = (
+            1j * 2 * np.pi / wavelength * (n_p - n_m)
+            * 4 / 3 * np.pi * radius**3
+        )
+
+        def scattered_field(mode, dx, refractive_index):
+            optics = Brightfield(
+                NA=1.2,
+                wavelength=wavelength,
+                resolution=dx,
+                magnification=1,
+                refractive_index_medium=n_m,
+                output_region=(0, 0, 256, 256),
+                padding=(128, 128, 128, 128),
+                return_field=True,
+            )
+            sphere = scatterers.MieSphere(
+                radius=radius,
+                refractive_index=refractive_index,
+                position=(128, 128),
+                z=0,
+                position_unit="pixel",
+                input_polarization=0,
+                output_polarization=0,
+                mode=mode,
+            )
+            return np.asarray(optics(sphere)())[..., 0]
+
+        for mode in ("hybrid", "geometric"):
+            for dx in (67e-9, 33.5e-9):
+                with warnings.catch_warnings():
+                    warnings.simplefilter("ignore")
+                    field = scattered_field(mode, dx, n_p)
+                    background = scattered_field(mode, dx, n_m)
+                ratio = (field - background).sum() * dx**2 / expected
+
+                # Sign: the phase delay of a denser sphere (both modes).
+                self.assertGreater(ratio.real, 0.5, msg=(mode, dx, ratio))
+                if mode == "hybrid":
+                    # Hybrid is exact in this limit: magnitude and phase.
+                    self.assertAlmostEqual(ratio.real, 1, delta=0.05,
+                                           msg=(mode, dx, ratio))
+                    self.assertAlmostEqual(ratio.imag, 0, delta=0.05,
+                                           msg=(mode, dx, ratio))
+
     def test_MieStratifiedSphere(self):
         scatterer = scatterers.MieStratifiedSphere(
             radius=(0.5e-6, 1.0e-6),