From fe2e0e803e57e11172dd7c4ac37817ba0f978bc6 Mon Sep 17 00:00:00 2001 From: edudc Date: Fri, 23 Oct 2026 12:51:18 +0200 Subject: [PATCH 1/4] fix sign and hybrid normalization --- deeptrack/optical/scatterers.py | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/deeptrack/optical/scatterers.py b/deeptrack/optical/scatterers.py index a7893f80f..7a6ca249c 100644 --- a/deeptrack/optical/scatterers.py +++ b/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( From 105375a9a4b1a941f1d1a8a6ea8b3d361a9cfac7 Mon Sep 17 00:00:00 2001 From: edudc Date: Fri, 23 Oct 2026 12:52:30 +0200 Subject: [PATCH 2/4] add test --- tests/test_scatterers.py | 52 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/tests/test_scatterers.py b/tests/test_scatterers.py index afc9fb784..bfff6a1f3 100644 --- a/tests/test_scatterers.py +++ b/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), From 95cd59cbe200c3e5d3e0fab4e16cb54fd83d7253 Mon Sep 17 00:00:00 2001 From: edudc Date: Fri, 23 Oct 2026 12:56:15 +0200 Subject: [PATCH 3/4] regenerated the gold values for test_dlcc test_3_B for the corrected Mie sign --- tests/test_dlcc.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/test_dlcc.py b/tests/test_dlcc.py index 85985d5db..63e2e9abb 100644 --- a/tests/test_dlcc.py +++ b/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(): From ec994f92ebeef04a8c61e8283dd1b61895cba993 Mon Sep 17 00:00:00 2001 From: edudc Date: Fri, 23 Oct 2026 17:10:59 +0200 Subject: [PATCH 4/4] remove plane wave from propagation. added test --- deeptrack/optical/optics.py | 4 ++-- tests/test_optics.py | 39 +++++++++++++++++++++++++++++++++++++ 2 files changed, 41 insertions(+), 2 deletions(-) diff --git a/deeptrack/optical/optics.py b/deeptrack/optical/optics.py index dc90b8ce0..a2040e3ed 100644 --- a/deeptrack/optical/optics.py +++ b/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( diff --git a/tests/test_optics.py b/tests/test_optics.py index ba5045ff5..6caf8e6d2 100644 --- a/tests/test_optics.py +++ b/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,